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Preface 


n 


This  work  covers  a  unique  combination  of  theory  and 
application  in  two  fields,  mathematics  and  computer  science.  It 
address  a  subject  with  a  wide  variety  of  uses  in  an  equally  wide 
variety  of  fields.  Thus  the  content  and  particularly  the  results  of 
this  work  should  appeal  to  a  large  audience. 

The  subject  I'm  referring  to  is  the  computation  of  four 
generalized  inverses.  The  significance  of  this  treatment  is  its 
improvement  over  existing  methods  for  computing  the  generalized 
inverses  of  a  matrix.  The  computational  method  developed  here  in 
theory  and  with  application  provides  a  simple,  direct,  and  unified 
approach  for  computing  four  generalized  inverse.  Probably  more 
importantly,  this  treatment  emphasizes  the  properties  of  each 
generalized  inverse  relative  to  its  computation  and  its  potential  use. 
For  some  readers,  hopefully  this  combination  will  provide  the 
missing  insight  needed  to  recognize  a  use  for  this  method  with  the 
proper  generalized  inverse  over  a  current  approach. 

Naturally,  with  an  effort  this  size,  there  are  a  number  of 
people  who  deserve  credit  and  a  word  of  thanks.  I  am  deeply 
indebted  to  my  thesis  advisor,  Dr  J.  Jones  Jr.,  for  his  interest, 
patience,  and  total  assistance  in  completing  this  effort.  Certainly  Dr. 
Jones  personifies  all  the  qualities  of  a  great  educator  and  scholar, 
especially  with  his  committment  and  enthusiasm.  1  also  thank  Capt 
Steve  Woffinden  for  his  comments  which  help  make  my  thesis  a 
more  understandable  document. 


ii 


Craig  F.  Murray 


Table  oL  Content? 


Page 


Preface . ii 

List  of  Figures .  v 

Abstact .  vi 


I.  Introduction .  1 

Background  .  1 

Characterizing  the  Generalized  Inverses  of  a  Matrix  ....  3 

Problem .  6 

Scope .  6 

Approach  and  Presentation .  7 

II.  Theory  of  Generalized  Inverses .  9 

Notation .  II 

Elementary  Row  and  Column  Operations . 12 

The  Generalized  Inverses  of  a  Matrix .  15 


O'v* 


III.  The  ST  Computation .  42 

The  ST  Method:  A  Synthetic  Approach .  43 

The  ST  Method  Adapted  to  Singular  Value  Decomposition.  65 

IV.  Application  to  the  Lyapunov  Matrix  Equation .  72 

The  Lyapunov  Matrix  Equation . 72 

The  Common  Solution . 82 

V.  Conclusions  and  Recommendations . 87 

Appendix  A:  Example  3 . 88 


iii 


Appendix  B:  Generalized  Inverse  Program .  96 

Appendix  C:  Graph  of  an  Aj  2  Versus  an  Aj  2  3^4 .  118 

Bibliography .  119 


m 


List  of  Figures 


Figure 

1.  The  Generalized  Inverses  of  a  Matrix . 

2.  Subset  of  Inverses . 

3.  Interchanging  Rows . 

4.  Multiplication  by  a  Scalar . 

5.  Multiplication  by  a  Scalar  and  Addition . 

6.  Matrix  Dimension  Table . 

7.  Dimensions  of  Augmented  Matrix  and  Identities 

8.  Partitioning  and  Dimensions  of  Submatrices  .  . 

9.  Five  Step  Summary . 

10.  Underdetermined  Case . 

11.  Algorithm  Comparison . 


V 


AFIT/G0R/MA/85D-2 


Abstract 

‘  " " 

This  w^k  examines  a  new  method  for  computing  four 
generalized  inverses  of  a  matrix.  This  method,  the  ST  method,  is 
based  on  the  careful  selection  of  a  sequence  of  matrix  multiplica¬ 
tions  and  partitionings  which  provide  a  new  foundation  for 
computing  four  generalized  inverses.  Central  to  this  approach  is  the 
partitioning  of  the  two  submatrices,  R  and  C,  where  the  product  of 
their  submatrices  will  give  the  generalized  inverses  of  interest. 

Thus  using  this  new  representation,  the  generalized  inverses  of  a 
matrix  can  be  computed  in  a  simple  and  direct  manner. 

In  this  work,  four  generalized  inverses  are  derived  and 
computed  in  a  systematic  manner  from  this  representation.  These 
results  are  strongly  tied  to  the  solution  of  the  matrix  equation 
Ax=b  where  the  general  solution  is  given  in  terms  of  this  new 
representation  and  computational  technique.  The  computational 
technique  is  presented  with  an  example  and  also  as  an  algorithm. 
Included  with  the  algorithm  is  an  analysis  of  its  computer 
implementation.  The  use  of  one  generalized  inverse  is  used  to  find 
the  solution  of  a  Lyapunov  matrix  equation. 


The  classical  inverse  of  a  matrix,  A,  is  another  matrix,  usually 
denoted  A“*,  which  satisfies  the  equation  AA'^  =  A~^A  =  I  .  For  the 
matrix  A  to  have  an  inverse,  in  this  sense,  one  important  restric¬ 
tion  must  be  met.  The  matrix  must  be  a  nonsingular  square 
matrix.  In  many  instances,  this  restriction  cannot  be  met.  For 
example,  in  systems  theory,  data  fitting,  statistics,  and  other  areas, 
rectangular  and  singular  matrices  arise  which  can  not  be  inverted 
using  the  classical  method  to  solve  the  respective  problem  (3;  5;  7:15; 
19).  Obviously,  in  these  cases  the  classical  inverse  can  not  be 
computed  so  what  must  be  done  to  solve  the  respective  problems? 

In  the  case  of  rectangular  or  singular  matrices  another  inverse,  or 
more  precisely  class  of  inverses  can  be  computed  and  used  to  find  a 
possibly  more  robust  solution.  This  class  of  inverses  is  called  the 
generalized  inverses  of  a  matrix. 

The  history  of  this  class  of  matrices  referred  to  as  generalized 
inverses  can  probably  be  succinctly  summarized  by  noting  five  key 
developments.  The  first  of  these  developments  is  creditted  to 
Fredholm  who,  it  is  believed,  introduced  the  concept  of  a  generalized 
inverse  in  1903  (3:4).  He  is  also  believed  to  have  introduced  the  term 
pseudoinverse  to  describe  his  new  concept  (3:4).  Following  Fredholm 


about  20  years  later  was  the  developments  set  forth  by  E.  H. 
Moore.  Moore  was  the  first  to  define  a  unique  generalized 
inverse  for  every  finite  matrix.  Moore  called  this  generalized 
inverse  the  general  reciprocal  (3:5).  While  the  progress  in  this  new 
field  continued  following  the  publication  of  Moore's  work,  the  next 
significant  development  came  another  30  years  later  when 
Bjerhammar  showed  the  relationship  of  a  generalized  inverse  to  a 
linear  system  of  equations  (3:5).  This  was  extended  four  years  later 
by  Penrose  who  showed  that  the  generalized  inverse,  defined  by 
Moore,  was  unique.  In  doing  this,  Penrose  defined  this  unique 
inverse  with  the  following  four  equations 


AXA  =  A 

(1.1) 

XAX  =  X 

(1.2) 

(AX)*  =  (AX) 

(1.3) 

(XA)*  =  (XA) 

(1.4) 

where  X  is  the  unique  generalized  inverse  of  the  matrix  A  (3:7). 

This  unique  generalized  inverse  is  popularly  referred  to  as  the 
Moore-Penrose  inverse.  After  Penrose’s  work,  the  next  major 
development  is  really  a  series  of  developments  starting  with  the 
work  of  Rao  who  defined  a  new  generalized  inverse  which  did  not 
satisfy  all  the  conditions  that  Penrose  introduced  with  Eqs  (1.1) 
through  (1.4)  (23:viii).  Instead  this  inverse  satisfied  a  combinations  of 
Penrose’s  equations.  Other  generalized  inverses  followed  which 
satisfied  fewer  or  different  combinations  of  the  four  equations  and 


were,  at  the  same  time,  more  robust.  The  major  draw  back  of 
these  developments  was  the  lack  of  a  direct,  unified,  and  numer¬ 
ically  stable  method  for  computing  the  generalized  inverses  whose 
existence  were  proved  theoretically.  This  problem  lead  to  the  fifth 
major  development  when  a  new  method,  the  ST  method,  was 
introduced  by  J.  Jones,  Jr.  in  1984.  This  latest  development  will  be 
documented  in  detail  in  this  work  emphasizing  its  computability. 

Characterizing  the  Generalized  Inverses  of  a  Matrix 


From  the  historical  notes  on  the  generalized  inverses  of 
matrices,  it  should  be  apparent  that  the  development  of  this  class  of 
matrices  suffered  a  rather  piecemeal  and  slow  growth.  This  has 
lead  to  an  equally  diverse  nomenclature  used  to  denote  these 
inverses.  Where  the  classical  inverse  of  the  matrix  A  was  simply 
denoted  as  a  A"^  and  defined  by  the  equation  AA"*  =  A‘^A  =  I,  the 
same  can  not  be  said  for  the  generalized  inverse  of  a  matrix. 

Instead  of  the  single  equation  found  in  the  classical  case,  there  is  a 
range  of  equations  for  the  generalized  inverses  of  a  matrix  where 
each  identifies  a  particular  property  a  given  inverse  satisfies.  This 
lead  to  the  wide  range  of  names  which  are  currently  used  to  refer 
to  these  inverses.  These  names  include  pseudoinverse, 

Moore-Penrose  inverse,  reflective  generalized  inverse,  g-inverse,  left 
weak  generalized  inverse,  and  right  weak  generalized  inverse. 
Additionally,  Ben-Israel  and  Greville  in  their  work  (3;18),  list  about 
24  other  names  used  to  describe  an  inverse  satisfying  some 


Figure  1.  The  Generalized  Inverses  of  a  Matrix 


combination  of  Penrose's  four  equations.  Thus  the  initial  problem  to 
overcome,  when  beginning  a  discussion  on  the  theory  of  generalized 
inverses,  is  choosing  a  suitable  method  for  describing  the  inverses  of 
concern.  One  means  of  overcoming  this  problem  is  by  adopting  a 
naming  convention  which  is  both  simple  and  direct  and  also 
incorporates  the  properties  of  a  given  inverse.  Starting  with  the 
given  matrix  A,  such  a  method  is  illustrated  in  figure  1  which 
graphically  depicts  a  partial  class  of  generalized  inverses.  This 
method  builds  from  the  given  matrix  A  up  to  the  inverse  satisfying 
all  four  of  Penrose's  equations,  Eqs  (1.1)  through  (1.4).  Thus  as  it 
travels  from  the  least  restricted  generalized  inverses  to  the  unique 
generalized  inverse,  Aj  2  3  ^ ,  the  combinations  which  arise  are 


4 


denoted  by  the  trailing  subscripts  where  each  subscript  refers  to 
one  of  the  four  equations  developed  by  Penrose.  The  subscripts 
listed  are  equated  to  the  conditions  satisfied  by  a  particular 
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Figure  2.  Subset  of  Inverses 

generalized  inverse.  To  prevent  confusion,  the  simple  relationship 
between  subscripts  and  the  four  equations  is  illustrated  below. 


Equation 
(AXA)  =  A 
(XAX)  =  X 
(AX)*  =  AX 
(XA)*  =  XA 


Corresponding  Subscript 
1 


w"-  W‘“  .V  wV  w>  1,-.  .V  -*  *.  .V  . . 


Again,  the  X  in  each  equation  denotes  the  appropriate  inverse. 
While  figure  1  shows  most  of  the  possible  generalized  inverses,  this 
work  will  concentrate  on  the  subset  of  these  inverses  shown  in 
figure  2  of  the  previous  page.  Again,  the  subscripts  following  each  A 
denote  which  of  Penrose's  four  equations  the  inverse  satisfies.  This 
nomenclature  will  be  used  throughout  this  work. 

Problem 

Current  methods  for  computing  the  generalized  inverses  of  a 
matrix  are  cumbersome,  and  many  are  also  numerically  unstable 
(19:247).  None  were  found  which  provide  a  unified  method  for 
calculating  more  than  one  generalized  inverse.  A  new  method  for 
computing  these  inverses,  the  ST  method,  is  a  direct  and  numer¬ 
ically  stable  technique  which  represents  a  marked  improvement 
over  previous  ones  in  use.  This  work  will  review  the  theoretical 
basis  of  the  ST  method,  and  go  on  to  define  both  the  characteristics 
of  its  computer  implementation  and  how  this  new  method  of 
computation  represents  an  improvement  over  past  methods. 

Scope 

In  this  work,  the  ST  method  is  considered  across  the  field  of 
complex  numbers.  For  the  treatment  of  theory,  this  means 
matrices  with  complex  elements  will  be  used.  In  terms  of  the 


technique's  computer  implementation,  this  translates  into  matrices 
of  constant  elements. 


In  describing  and  defining  the  technique’s  computer  implemen¬ 
tation,  three  areas  will  be  examined.  First  the  numerical  properties 
of  this  technique  will  be  treated.  This  will  be  followed  by  a  time 
and  space  analysis  of  the  algorithm  arising  from  the  ST  method. 
Finally,  this  technique  will  be  compared  to  other  current  techniques 
to  establish  its  claimed  improvements. 

Approach  and  Presentation 

Chapter  11  begins  the  review  of  the  theoretical  basis  of  the  ST 
method  for  computing  the  generalized  inverses  of  a  matrix  with 
complex  elements.  The  theory  begins  with  a  new  computation 
representation  and  continues  to  present  a  systematic  development 
of  four  generalized  inverses.  Each  theory  is  developed  with  as  much 
detail  as  possible  to  insure  maximum  clarity.  An  example  is  given 
early  during  this  treatment  to  motivate  and  assist  the  understand¬ 
ing  of  the  theorems  presented. 

Chapter  111  applies  the  theorems  developed  in  chapter  11  to  the 
realm  of  the  computer.  Here  a  detailed  description  of  the  technique 
is  given  in  an  algorithmic  manner  versus  the  theoretical  basis  of 
chapter  II.  In  this  regard,  the  algorithm  for  the  ST  method  is 
examined  to  determine  its  numerical  characteristics  and  analyze  its 
computer  time  and  space  requirements.  This  is  followed  by  an 
introduction  to  other  techniques  for  computing  the  generalized 


inverses  of  a  matrix,  and  these  are  compared  to  the  ST  method  in 
terms  of  numerical  stability  and  computer  requirements.  Finally, 
the  ST  method  is  adapted  to  an  existing  algorithm  available  in 
many  computer  math  packages  to  compute  an  Aj  2  3  4  generalized 
inverse. 

Chapter  IV  concludes  the  heart  of  this  work  by  extending  the 
comparisons  developed  in  chapter  III  to  an  actual  application  of  the 
ST  method.  In  this  regard,  the  ST  method  is  used  to  find  the 
solution  of  a  Lyapunov  matrix  equation.  To  do  this,  additional 
theorems  are  developed  to  show  how  the  ST  method  is  applied  to 
the  Lyapunov  matrix  equation.  Finally,  an  Al,2  generalized  inverse 
is  developed  in  the  common  solution  for  a  Lyapunov  matrix 
equation. 

Chapter  V  completes  this  work  with  a  summary  of  its 
important  points  and  recommendations  for  further  study  in  this 


II  Theory  of 


Inverses 


As  mentioned  earlier,  the  theory  of  generalized  inverses  has 
been  introduced  under  a  variety  of  names,  including  the  pseudo¬ 
inverse  and  the  Moore-Penrose  inverse.  Common  to  each  approach, 
no  matter  what  the  name,  is  an  emphasis  on  the  generalized 
inverse  which  satisfies  the  four  equations; 


AXA  =  A 

(2.1) 

XAX  =  X 

(2.2) 

(AX)*  =  AX 

(2.3) 

(XA)*  =  XA 

(2.4) 

Subsequent  to  such  an  introduction,  this  common  approach  con¬ 
tinues  by  introducing  a  variety  of  properties  arising  from  the 
pseudoinverse  or  Moore-Penrose  inverse  satisfying  Eqs  (2.1)  through 
(2.4),  and  then  gradually  moves  toward  a  more  complete  discussion 
of  the  other  generalized  inverses  of  a  matrix.  Such  an  approach 
concludes  by  finally  presenting  a  method  for  computing  the 
generalized  inverses  again,  emphasizing  the  one  inverse  satisfying 
Eqs  (2.1)  through  (2.4). 

The  approach  to  the  theory  of  generalized  inverses  presented  in 
this  work  will  be  different.  Instead  of  the  common  path  taken  by 
many  authors  in  the  past,  this  approach  presents  a  systematic 


development  of  four  generalized  inverses  with  their  respective 
properties  and  their  computability.  Additionally,  this  discussion  of 
the  generalized  inverses  of  a  matrix  is  strongly  oriented  toward 
solving  the  system  of  equations  often  denoted  as  Ax=b  where  the 
matrix  A  does  not  have  a  classical  inverse.  Central  to  this 
orientation  is  an  emphasis  on  computation  which  will  be  presented 
hand  in  hand  with  theory.  This  is  intended  to  provide  both  an 
understanding  of  theory  and  the  application  of  the  theory.  Thus 
when  the  equivalent  of  the  psuedo-  inverse  or  Moore-Penrose 
inverse  is  finally  reached,  its  properties  and  computability  arises 
more  directly.  Unfortunately  this  approach  shares  one  common 
problem  with  its  predecessors,  the  problem  of  nomenclature.  Since 
the  means  of  describing  a  generalized  inverse  of  a  matrix  varies,  as 
noted  in  chapter  1,  a  few  additional  comments  are  devoted  to  more 
completely  defining  the  notation  used  in  this  work.  Other  remarks 
regarding  notation  used  throughout  this  work  are  also  included 
at  this  point  to  prevent  confusion  and  emphasize  the  important 
characteristics  of  the  material  which  follows.  Additionally,  since  this 
work  is  geared  to  the  practical  adaptation  of  a  generalized  inverse 
to  solving  real  world  problems,  a  section  will  be  devoted  to  the  basic 
operations  used  to  compute  the  generalized  inverses,  namely 
elementary  row  and  column  operations.  For  the  reader  unfamiliar 
with  these  operations,  this  review  will  be  helpful  for  understanding 
the  computations  involved  in  finding  the  generalized  inverses  of 
concern  here,  For  the  reader  familiar  with  row  and  column 
operations,  their  mention  here  should  provide  an  insight  to  the 
directness  of  the  computation  to  be  discussed. 


Notation 


The  remarks  concerning  the  notation  used  throughout  this 
work  deal  with  both  the  exact  expressions  used  to  denote  the 
generalized  inverses  of  interest  and  the  mathematical  expressions 
used  to  derive  these  inverses.  Both  are  discussed  below. 

As  mentioned  in  chapter  1,  the  generalized  inverses  of  a  matrix 
have  been  denoted  in  a  variety  of  ways.  In  this  work  the  following 
representation  will  be  used 

^ij.k,l 

where  the  subscripts  denote  the  equations  the  given  generalized 
inverse  satisfies.  Thus  a  generalized  inverse  satisfies  Eqs  (2.1)  and 
(2.2)  would  be  denoted  as 

^1,2 

A  generalized  inverse  satisfying  all  four  equations  would  be  denoted 
as 

^1,2, 3.4 

This  notation  provides  an  easy  and  direct  method  for  classifying  and 

describing  the  generalized  inverses  of  a  given  matrix. 

The  second  notational  concern  deals  with  the  mathematical 
expressions  used  throughout  this  work.  First,  the  symbol,  Q,  is  used 
to  denote  the  field  of  complex  numbers,  and  denotes  the 

vector  space  of  m  x  n  complex  matrices  over  Q.  Second,  matrices 
are  represented  by  capital  letters,  and  the  small  letters  b,  w,  x,  y, 
and  z  denote  vectors  in  this  work.  Thus,  A,  would  denote  a  matrix 
A  of  unspecified  dimensions.  Third,  an  asterisk  denotes  the  conjugate 


transpose  of  a  matrix.  Thus  the  conjugate  transpose  of  A  would  be 
denoted  as  A*.  The  final  comment  on  notation  again  involves  the 
representation  of  matrices.  When  depicting  a  large  matrix  with 
multiple  submatrices,  the  matrix  will  be  enclosed  in  brackets,  and 
the  submatrices  denoted  as  capital  letters.  In  some  cases,  the 
submatrices  will  also  be  separated  by  partitioning  lines.  This 
practice  is  only  used  in  a  few  cases  to  emphasize  the  composition  of 
a  matrix.  Normally,  the  partitioning  of  a  large  matrix  will  be 
denoted  using  the  submatrices  depicted  in  capital  letters,  sufficient 
spacing,  and  the  enclosing  brackets.  When  numeric  values  are  listed 
as  the  components  of  a  matrix  the  meaning  of  the  above 
conventions  changes  slightly.  If  the  value  is  a  zero,  then  it  may 
represent  the  real  number  zero  or  the  zero  matrix.  While  this  may 
sound  confusing,  the  context  of  such  an  occurrence  will  clarify  its 
meaning.  Other  occurrences  of  numeric  values  in  matrices  do  depict 
individual  matrix  elements. 

Elementary  Row  and  Column  Operations 

In  the  next  section,  theorem  2-1  will  provide  the  basis  for 
computing  the  generalized  inverses  of  a  matrix.  To  understand  how 
to  utilize  this  result,  in  a  computational  sense,  a  knowledge  of 
elementary  row  and  column  operations  is  required.  Further  an 
understanding  of  these  operations  will  help  develop  an  appreciation 
for  the  computationally  direct  and  simple  process  used  to  obtain 
these  generalized  inverses.  Thus  as  either  a  review  or  quick 


Figure  3.  Interchanging  Rows 


overview,  the  following  summary  of  elementary  row  and  column 
operations  will  help,  at  a  minimum,  appreciate  the  computational 
aspects  of  the  processes  at  the  heart  of  this  work. 

Interchanging  Rows  ql  Columns.  Given  a  matrix,  any  complete 
row  of  elements  may  be  interchanged  with  the  corresponding 
elements  of  any  other  complete  row  of  the  same  matrix.  The  new 
matrix  is  equivalent  to  the  original  matrix.  Using  an  arbitrary 
matrix,  this  operation  is  demonstrated  in  figure  3  where  the  first 
and  third  rows  are  interchanged.  The  same  conceptual  operation 
can  also  be  performed  on  columns.  Now  instead  of  interchanging  the 
corresponding  elements  of  two  rows,  the  corresponding  elements  of 
two  columns  are  interchanged.  Had  this  operation  been  applied  to 
the  matrix  in  figure  3,  interchanging  columns  one  and  three  would 
not  have  produced  a  matrix  where  the  first  column  consisted  of  a 
three,  two,  and  one.  Interchanging  columns  one  and  two  would  not 


have  made  a  noticeable  difference  since  each  column  is  identical. 


Multiplication  iif  a  Row  ql  Column.  Given  a  matrix,  the 
elements  of  any  row  in  the  matrix  can  be  multiplied  by  a  scalar. 
The  resulting  matrix  is  equivalent  to  the  original  matrix.  Again, 
using  the  same  arbitrary  matrix,  this  operation  is  demonstrated 
below  in  figure  4  where  row  1  is  multiplied  by  scalar  k=3.  This 
same  operation  may  be  performed  on  columns  as  well  as  rows. 

Thus  had  the  elements  of  column  1  been  multiplied  by  k=3  the 
resulting  column  would  contained  3,  6,  and  9  from  top  to  bottom 
respectively. 

Multiply  a  Row  ql  Column  and  Add-  Using  the  same  matrix 
from  the  previous  illusrtation,  a  row  is  again  multiplied  by  a  scalar 
and  the  resulting  new  row  is  now  added  to  one  of  the  remaining 
rows  of  the  matrix.  As  before,  the  addition  is  performed 


\V, 


1  1  I 

2  2  2 
3  3  3 

Original  Matrix 

1  1  l' 

2  2  2 
0  0  0 

Matrix  After  Multiplying  Row  1  by  Scalar  k=(-3) 
and  then  Adding  Row  1  to  Row  3 


Figure  5.  Multiplication  by  a  Scalar  and  Addition 


on  only  corresponding  elements  of  each  row.  The  results  of 
multiplying  row  1  by  k=(-3)  and  then  adding  this  new  row  to  row 
3  are  shown  in  figure  5.  Note,  row  1  is  unchanged  by  this  operation 
and  row  3  is  eliminated.  The  same  operation  can  be  performed  on  a 
column  basis.  Using  column  operation  and  the  same  original  matrix 
from  figure  5,  column  1  could  be  used  to  eliminate  the  elements  of 
column  3  by  multiplying  its  elements  by  k=(-l)  and  then  adding 
corresponding  elements  of  columns  1  and  3.  Had  this  same  operation 
been  performed  on  the  resepective  columns  of  figure  3,  only  the 
second  element  in  column  three  would  have  been,  eliminated. 


Ih£  Generalized  Inverses  qL  a  Matrix  (10) 


The  systematic  development  of  the  four  generalized  inverses 
begins  with  the  construction  of  a  new  and  equivalent  represen- 
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tation  of  the  arbitrary  matrix  A  which  has  dimensions  m  rows  and 
n  columns.  The  starting  point  of  this  systematic  development  is 
then  the  matrix  A  augmented  with  two  identity  matrices  of 
appropriate  dimensions  and  shown  in  the  following  equation 


Using  this  representation,  the  the  four  generalized  inverses  of 
interest  will  be  developed  in  theory  and  in  practice.  This 
representation  and  consequently  the  generalized  inverses  derived 
from  it  will  also  be  closely  tied  to  the  general  solution  of  a  system 
of  equations.  This  is  all  accomplished  through  seven  theorems  and 
one  example.  Theorem  2-1  begins  by  developing  the  new 
representation  from  Eq  (2.5).  This  is  immediately  followed  by  an 
example  which  demonstrates  how  each  of  the  four  generalized 
inverses  of  interest  is  computed  from  this  new  computational 
framework.  Theorem  2-2  extends  this  example  by  introducing  how 
the  general  solution  of  a  system  of  equations  also  arises  from  this 
computational  framework  in  terms  of  the  generalized  inverses. 
Theorems  2-3,  2-4,  2-5,  and  2-7  develop  the  existence  of  each  of  the 
generalized  inverses  computed  from  this  new  representation.  Finally 
theorem  2-6  develops  the  general  solution  of  a  system  of  equations 
in  a  manner  similar  to  theorem  2-2  but  where  the  submatrix  M  is 
not  defined  based  on  the  rank  of  A. 

Theorem  2d  For  any  given  matrix  A  e  g’^xn^  there  exist  two 


■  ^  ^  j|  ' 

m^.  JkB  aA  M-t. 
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nonsingular  matrices,  R  e  q  ^  gnxn^ 

following  pair  of  equivalent  matrices  can  be  constructed 


A  I 
I  0 


and 


'k 

0 

t' 

0 

0 

IVI 

N 

0 

(2.6) 


Proof.  For  any  matrix  A  e  there  exist  nonsingular 

matrices  R  e  gf^^xm  q  ^  gJ^xn 


Ir  0 
0  0 


(2.7) 


where  1^.  is  an  identity  matrix  of  dimension  equal  to  the  rank  of 

matrix  A.  Eq  (2.7)  is  a  well  established  result,  and  this  result  is  used 
as  It  arises  in  the  following  matrix  multiplications 


'p.  o' 

A  r 

C  ') 

RA  r' 

C  0 

~ 

F:AC  P: 

I. 

I  0. 

.0  I, 

I  0. 

0  I 

C  0 

Now  matrices  R  and  C  are  now  defined  as  follows 


R 


::ind  L 


[S  N] 


(2.y,i 


Using  these  new  forms  of  R  and  C,  the  desired  result  is  obtained 
from  Eq  (2  8)  giving  the  following 


RAC  R 
r.  0 ! 


1,-1  0  T 
0  10  M 
SIN  0 


(2.10) 


At  this  point,  the  definition  of  matrices  R  and  C  in  Eq  (2.8)  may 
seem  somewhat  arbitrary,  but  these  new  forms  will  play  an 
important  role  in  computing  the  various  generalized  inverses  of  the 
given  matrix  A  c  This  representation,  developed  in  theorem 

2-1,  will  allow  the  partitioning  of  A  into  the  matrices  S,  T,  M,  and  N 
using  a  scheme  reminiscent  of  the  classical  approach  to  matrix 
inversion.  These  new  matrices  are  then  the  foundation  for  comput¬ 
ing  the  four  generalized  inverses  of  A  as  discussed  above.  The 
following  example  will  help  illustrate  how  the  S,  T  M,  and  N 
matrices  are  obtained,  and  also  how  these  matrices  will  form  the 
foundations  for  computing  the  generalized  inverses  of  interest. 

Example  1.  Given  the  matrix  A  where 


1.000  2.000  3.000  4.000 
5  0  0  0  1 0 . 0  0 0  15. 0 0 0  2 0  0 0 0 


(2.11) 


Augment  A  with  an  identity  matrix,  beside  and  below  it,  to  give 
the  representation  of  Eq  (2.12); 


C'l. 

1  ’  I'i  ’ 

L  ‘  '-'j 


1  1  I 

'  1 

l.i 

4  0 

C  I'i  1  0 

0  15 

0 

20  0 

1  ' '  0 

'•  ■’  () 

(i 

0  (' 

I'I  (i  1 

'  0 

!  1 

0  0 

1,1  0  I) 

0  1 

(; 

0  0 

1'  1  ■'<  !'■» 

11  ri 

11 

I  1  f  I  u  n 
)  0  0  1  0 


(2  12) 


Now  using  elementary  row  and  column  operations,  reduce  A  to  the 
following  form  where  I  is  an  identity  matrix  of  dimension  equal  to 
the  rank  of  A  and  the  identity  matrices  of  Eq  (2.5)  are  not  shown; 

I  0 
0  0 

From  this  point  on,  the  identity  matrix  1,  arising  from  the  reduction 
of  A,  will  be  referred  to  as  Ij.  as  introduced  in  theorem  2-1.  The 

form  of  Eq  (2.12)  should  look  familiar.  It  was  the  starting  point  of 
theorem  2-1,  and  it  was  eventually  transformed  it  to  an  equivalent 
representation  which  exhibited  the  submatrices  S,  T,  M,  and  N. 

This  representation,  Eq  (2.10),  will  be  obtained  in  this  example  by 
using  elementary  row  and  column  operations  to  reduce  the  matrix 
A  to  the  identity  matrix  1^..  As  A  is  reduced,  the  desired  submatric¬ 
es  will  arise  on  the  two  augmented  identity  matrices  similar  to  the 
classical  method,  where  the  original  matrix  is  augmented  with  an 
identity  matrix  and  then  original  matrix  is  reduced  to  an  identity 
matrix.  In  the  classical  case,  the  result  of  the  row  operations  on  the 
identity  matrix  produce  the  classical  inverse,  but  here  both  row 
and  column  operations  will  produce  four  new  submatrices,  two  of 
which  will  be  used  to  compute  the  generalized  inverse  of  interest. 

The  process  begins  by  first  performing  the  desired  row 
operations  to  eliminate  the  second  row  of  A  as  illustrated  in  Eq 
(2.13).  As  the  second  row  of  A  is  eliminated,  the  results  of  this  row 
operation  begin  to  define  the  submatrices  of  interest.  Next, 


m 


1.0  2.0  3.0  4.0  0.0  0.2 

0.0  0.0  0.0  0.0  1.0 -0.2 

1.0  0.0  0.0  0.0 

0.0  1.0  0.0  0.0 

0.0  0.0  1.0  0.0 

0.0  0.0  0.0  1.0 


(2.13) 


performing  the  required  column  operations,  A  is  reduced  to  as 
shown  in  Eq  (2.14); 


1.0  0.0  0.0  0.0  0.0  0.2 

0.0  0.0  0.0  0.0  1.0 -0.2 

1 .0  -2.0  -3.0  -4.0 

0.0  1 .0  0.0  0.0 

0.0  0.0  1.0  0.0 

0.0  0.0  0.0  1 .0 


(2  14) 


The  augmented  version  of  A,  Eq  (2.12),  has  now  been  transformed  to 
the  equivalent  representation  exhibiting  the  submatrices  S,  T,  M, 
and  N.  These  are  illustrated  more  explicitly  in  Eq  (2.15): 


1.010.0  0.0  0.0 !  0,0  0.2 


0  0  i  0.0  0.0  0.0 


2.0  -3.0  -4.0 
10  0.0  0.0 
0.0  1 .0  0.0 
0  0  0.0  1  0 


I  .ouo 


I 


-2.000-3.000-4.000 
1 .000  0.000  0,000 
0.000  1 .000  0.000 
0.000  0.000  1 .000 

T  =  [0  000  0.200] 

M  =  [1.000-0.200] 

With  the  S,  T,  M,  and  N  matrices  computed,  the  next  objective 
of  this  example  is  to  form  the  generalized  inverses  of  A  using  these 
new  submatrices.  To  do  this  one  simple  rule  is  required.  The 
generalized  inverse  is  equal  to  the  product  of  the  submatrices  S  and 
T.  Using  this  rule  and  the  notation  introduced  earlier,  the  Aj  2 

generalized  inverses  will  be  computed  directly  from  the  result  in  Eq 
(2.15). 


The  Aj  2  generalized  inverse  is  formed  from  the  product  of  the 
kV,  matrices  S  and  T.  Thus,  the  multiplication  in  Eq  (2.16)  directly  and 

simply  obtains  the  Aj  2  generalized  inverse  of  A; 


■■r. 

<• 


■>. 


'  1 .000' 

0  000  0.200 

0.000 

[0.000  0.200]  = 

0.000  0.000 

1,1  i,_(i.jg 

(.1  Cl 00  0 

0  000 

0.000  0  000 

where  Aj  2  satisfies  the  first  two  conditions  defined  by  Penrose; 
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I  wmVji  ^  n  ^  *  ta.  Ji  .  /I  ^  1  s  \ 


A(ST)A  =  A 


(2.17) 


(ST)A(ST)  =  ST 


(2.18) 


While  the  Aj  2  generalized  inverse  was  computed  directly  from 
Eq  (2.15),  obtaining  the  Aj2  3  and  k^2A  generalized  inverses,  as  one 

might  expect,  requires  a  little  more  work.  Both  inverses  are  still  the 
product  of  S  and  T,  but  now  the  form  of  S  and  T  change  depending 
on  which  inverse  is  desired.  For  A|  2  3  each  row  of  T  is 

orthogonalized  to  the  rows  of  M.  For  the  Aj  2  4  generalized  inverse, 

each  column  of  S  is  orthogonalized  to  the  columns  of  N.  When  both 
orthogonalization  processes  are  performed,  S  and  T  appear  as  shown 
below  (To  orthogonalize  the  rows  or  columns,  the  Gram-Schmidt 


1 .000  0.000 

0.000  0.000  0. 

038  0.192' 

0.000  0.000 

0.000  0.000  1. 

000-0.200 

0.033-0.077- 

0. 1 76  -4.000 

0.067  1.000 
0. 1 00-0.23 1 

0.0 Oo  u.OOu 

1.000  0.000 

0.133-0.308- 

0.706  1.000 

orthogonalization  process,  was  used.  A  more  detailed  explanation  of 
this  process  is  offered  in  chapter  3). 

Computing  the  Aj  2  3  generalized  inverse  requires  a  step  back 
from  the  fully  orthogonalized  representation  shown  above  to  Eq 
(2.15).  To  computed  the  Aj2  3  inverse  requires  the  orthogonalized  T 
submatrix  and  the  original  S  submatrix  from  Eq  (2.15).  When 


i 


^^^ /* 


performed,  the  multiplication  appears  as  illustrated  in  Eq  (2.19): 


'1  000] 

0.038  0.192 

0.000  ^ 

0.000  0.000 

0.000 

o.oooj 

0.000  0.000 

0.000  0.000 

(2.19) 


where  S  is  the  original  S  submatrix  from  Eq  (2.15)  and  T  is  now  an 
orthogonalized  version  of  the  T  from  Eq  (2.15).  This  Aj2  3  generalized 

inverse  satisfies  Eqs  (2.17)  and  (2.18)  plus  the  additional  Penrose 
condition  detailed  in  Eq  (2,20): 


(ST)(A)»  =  (ST)A 


(2.20) 


The  process  for  computing  the  Aj  2  4  generalized  inverse 

proceeds  in  a  similar  manner.  Instead  of  the  orthogonalized  version 
the  the  T  submatrix  and  the  original  S  submatrix,  the  opposite  of 
each  is  used.  The  resulting  multiplication  follows  in  Eq  (2.21) 


(2.2(;?) 


'0.033' 

0.000  0.007' 

2,4=  •3T  = 

0.067 
0. 100 

[0  0  0.21  = 

0.000  0.013 
0.000  0.020 

0  1  ^ 

A  AAA  A  A’"!*"? 

U.UUv  U.U-L  / 

Again,  along  with  Eq  (2.17)  and  (2,18),  the  Aj2  4  generalized  inverse 
satisfies  Eq  (2.22): 

(ST)*(A)  =  (ST)A  (2.22) 


(2.22) 


With  three  out  of  a  possible  four  choices  used,  there  is  only  one 
logical  combination  of  the  S  and  T  submatrices  left  for  computing 
the  final  generalized  inverse  of  interest,  the  Aj  2  3  4  inverse. 

The  2  3  4  verse  is  obtained  in  a  similar  manner  using  the  S 

and  T  submatrices,  but  now  both  orthogonalized  versions  are 
multiplied  together.  In  this  example  the  result  is  obtained  by 


0  033 

4.2, 3, 4=  ST=  10.038  0.192] 

u.luu 

0,133 


0  001  0  006 
0.003  0  013 
0  004  0.019 
0.005  0  026 


(2.2o) 


The  Aj  2  3  4  inverse  satisfies  all  the  conditions  satisfied  by  the  lower 
inverses  plus  one  last  condition  making  it  the  unique  inverse: 


(A(ST))»  =  A(ST) 


(2.24) 


A  second  more  detailed  example  is  provided  in  appendix  A.  This 
example  was  taken  from  pages  340  and  341  of  Noble's  Applied 
Linear  Algerbra  text,  reference  21.  All  four  generalized  inverses  are 
computed  and  each  of  the  four  conditions  are  tested  to  verify  the 
results. 

Theorem  2-1  defined  a  new  starting  point  for  computing  the 
generalized  inverses  of  a  matrix.  In  this  respect,  theorem  2-1 
provided  both  a  new  method  of  representing  a  matrix  and  a  new 
means  of  partitioning  a  matrix  to  a  form  well  suited  for  computing 
the  generalized  inverses  of  interest.  Using  this  result,  example  1  took 


a  2x4  matrix  and  reduced  it  to  the  form  developed  in  theorem  2-1, 
the  new  starting  point  for  computing  the  generalized  inverses  of 
concern.  The  process  of  reducing  the  matrix  to  the  desired  form 
was  computationally  simple  and  direct  using  only  elementary  row 
and  column  operations.  From  this  starting  point,  the  process  was 
expanded  to  compute  four  generalized  inverses  where  each  inverse 
satisfies  one  or  more  of  the  four  conditions.  Combined,  theorem  2-1 
and  example  1  provide  a  starting  point  and  insight  into  how  the 
result  and  process  can  be  put  to  further  use.  Theorem  2-2  will 
begin  to  expand  on  this  by  extending  the  new  representation  and 
computational  technique  to  the  framework  of  the  equation  Ax=b  . 

Theorem  2-2.  The  system  of  equations  defined  by 

Ax  =  b 
where 

A  e 

X 

b 

has  a  solution  x  if  and  only  if  Mb  =  0  and  the  general  solution  of 
Ax=b  can  be  given  by 

^general  solution  "  (2.25) 

where  S,  T,  M,  and  N  are  the  matrices  defined  in  theorem  2-1  and 
z  is  an  arbitrary  variable  of  appropriate  dimensions. 


Proof.  In  theorem  2-1,  matrices  R  and  C  were  partitioned  as 


R  =  J  and  C  =  [S  N] 
Ivl 


(2.9) 


Additionally  R  and  C  are  nonsingular  as  defined  in  theorem  2-1.  By 
using  this  partitioning  along  with  matrix  multiplication,  theorem  2-1 
provided  the  following  relationship 


RAC  R 
C  0 


If  I  0  T 
^  IVI 


(2,10) 


From  these  matrices  and  their  submatrices,  the  following  chain  of 
implications  arise 

Ax  =  b  has  a  solution  x, 

if  and  only  if  RAx  =  Rb  has  a  solution  x, 

if  and  only  if  RACy  =  Rb  has  a  solution  y  and  x  =  Cy, 

if  and  only  if  b  has  a  solution  v  = 

.0  0j[zj  [mJ  [z. 


and  A  =  [S  N] 


if  and  only  if 


'ir  W' 

Tb' 

has  a  solution  y  = 

'v/' 

0 

.Mb 

»7 

and  X  =  [SW  +  Nz]  , 

if  and  only  if  1^,  W  =  Tb,  0  =  Mb,  and  x  =  [SW+Nz] 

Thus  from  the  above  chain  of  implications,  Mb=0  is  a  consistency 
condition  and  determines  if  the  the  equation 


X  =  SW+Nz  (2.26) 

has  a  solution.  Further,  Eq  (2.26)  can  be  rewritten  in  the  following 
form  to  capitalize  on  the  results  of  theorem  2-1,  and  obtain  the 
results  stated  in  theorem  2-2: 

X  =  S(Tb)+Nz  (2.27) 


Corollary  1.  By  choosing  b=0  in  theorem  2-1,  the  columns  of  N 
form  a  basis  for  the  null  space,  TifA),  where  Ax=b  has  the  general 
solution 

^general  solution  “  ^  (2.28) 


Again,  z  must  be  of  appropriate  dimensions  and  M*b=0  is  a 
consistency  condition 


Proof.  By  choosing  ,  the  columns  of  N  span  'n(A),  Also  the 
columns  of  A  come  from  the  nonsingular  matrix  C.  Therefor  the 
columns  of  N  are  linearly  independent  and  hence  form  a  basis  for 


With  a  growing  understanding  of  the  ST  method  from  theorem 
2-1  and  2-2,  nothing  yet  has  been  said  about  how  the  product  of 
the  S  and  T  submatrices  satisfies  combinations  of  the  four  equations 
defined  by  Penrose,  The  next  three  theorems  plus  theorem  2-7  will 
prove  the  existence  of  each  generalized  inverse  of  concern  in  this 
work  relative  to  the  results  of  theorem  2-1  and  theorem  2-2. 


Theorem  2-3.  From  the  representation  developed  in  theorem 
2-1  and  the  process  described  in  example  1,  the  matrix  ST  is  a  Aj 

generalized  inverse  of  A  where  Aj  e  and  A  €  Q^^^n 

generalized  inverse  Aj  satisfies  the  following  conditions: 


(ST)A(ST)  =  ST  =  Aj 


Proof.  Using  the  R  and  C  matrices  defined  earlier  in  Eq  (2.7), 
the  matrix  multiplication  of  RAC  can  be  written  slightly  differently 
than  its  form  in  Eq  (2.8).  The  RAC  multiplication  can  now  be 
extended  with  the  four  submatrices  S,  T,  M,  and  N  providing 
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Figure  6.  Matrix  Dimension  Table 


The  next  step  is  to  begin  a  long  chain  of  multiplications  and 
substitutions  which  can  be  deceiving  when  using  matrices  without 
clearly  depicted  dimensions.  Figure  6  depicts  the  correct  dimension 
of  each  matrix  appearing  in  the  subsequent  multiplications. 
Reference  figure  6  if  any  confusion  arises  regarding  the  conformity 
of  a  multiplication  or  substitution.  Now  applying  the  multiplications 
introduced  in  Eq  (2.29) 


[S  N] 


TA 

MA 


[s  m 


(2.30) 


Ti  1 

1  [^1  M]  =  ! 

TAc^  TAM' 

llVlAj  1 

MAS  MAN 

(2  31) 
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(ST)A(ST)  =  ST  =  Aj 

holds,  the  result  from  Eq  (2.31)  are  substituted  back  into  Eq  (2.29) 
giving 

'ir  ol  [tAS  TAnI 

0  oj  [ivlAlS  man 

This  relationship  implies 

TAS  = 

Thus 

(ST)A(ST)  =  S(TAS)T  =  S(Ir)T  =  ST  =  Aj 

■ 

While  the  Aj  generalized  inverse  was  not  computed  in  example 
1,  it  will  be  used  here  as  an  intermediate  step  to  reaching  the  Aj  2 

generalized  inverse.  Thus  the  next  step  will  be  to  show  ST  also 
satisfies  the  second  condition  introduced  in  Eq  (2.2).  Once  this  is 
proved,  ST,  in  this  case,  will  be  an  Aj  ^  generalized  inverse  by 
definition  since  it  satisfies  Penrose’s  first  two  conditions. 


2-4.  The  ST  matrix  obtained  in  theorem  2-3  is  also  an 


A2  generalized  inverse  satisfying 


A(ST)A  =  A 


Proof.  Recalling  from  above 


PiL. 


Ir-  0 
0  0 


Since  R  and  C  are  nonsingular,  it  follows  that 


R 


:.-1  [Ir 
[0 


0 


Using  the  new  definition 


Ir- 


Eq  (2.33)  can  be  expressed  as 


Ir  0 
0  0 


(2.17) 


Now  consider  the  following  using  the  previously  defined  properties  of 


R  =  [3mp 


(i.otO 


If  0 

t' 

0  0 

M 

==  [Sir  0] 


=  S(Ij,)T  =  ST 

Using  t’  ese  and  previous  results  to  substitute  for  A  and  ST  gives 

A(ST)A  =  (R-1bC'1)(CB*R)(R-‘BC‘*) 

=  R-‘B(C'1C)B*(RR"^)BC-1 
=  R'^BB*BC'^ 

=  R'^C'l 
=  A 

Finally  combing  this  result  on  top  of  the  previous  result,  the  proof 
of  theorem  2-A  follows  with 

(ST)A(AT)  =  (CB*R)(R-^BC-1)(CB*R) 

=  CB*(RR-‘)B(C'‘C)B*R 
=  CB*R 
=  ST 
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Once  again,  since  the  same  ST  matrix  was  used  to  obtain  the 
generalized  inverses  and  A2  in  theorem  2-3  and  theorem  2-4 
respectively,  ST  satisfies  both  conditions  and  is  then  by  definition  an 
Aj  2  generalized  inverse.  The  same  type  of  argument  will  be  used  to 

prove  the  existence  of  the  Aj  2  3  4  generalized  inverse.  This  begins 

with  theorem  2-5  below  and  concludes  with  theorem  2-7.  Theorem 
2-6  interrupts  the  process  to  provide  the  additional  results  needed 
to  prove  the  existence  of  an  Aj  2  4  generalized  inverse. 

Theorem  2-5.  If  MT*=0,  then  the  matrix  ST  is  an  Aj  2  3 
generalized  inverse  of  the  matrix  A. 

Proof.  It  must  be  shown  that  the  following  equations  hold  for 
any  given  A  €  and  Aj  2^3  e 


^  ^1,2,3  ^  ^ 

(2.37) 

^1,2,3  A  ^1,2,3  "  ^1,2,3 

(2.38) 

^^^1,2,3^  "^^^1,2,3^ 

(2.39) 

Let  A'  =  ST.  It  is  sufficient  to  show  that  AA~  is  symmetric  since  by 
the  results  of  theorem  2-4  it  is  already  know  that 

A  A'A  =  A 
A'A  A'  =  A" 
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v.v.v 


Then  in  order  to  show  any  matrix  is  symmetric,  it  is  sufficient 


show  that 


(A  AKA  A’)*  =  (A  A") 


(2.40) 


In  other  words,  A  is  symmetric  if  and  only  if  A  =  A  .  Now  to  show 
that  (AA”)  is  symmetric,  begin  by  showing  (AA')(AA')*  is  symmetric. 
First  note  that  (PQ)*=  Q*P*  holds  for  matrix  multiplication  of  P  and 
Q.  Then  taking  the  transpose  of  (AA")(AA~)*  as  with  P  and  Q  gives 


[(AA‘)(AA')*]*  =  ((AA')*)*(AA‘)*  =  (AAKAAl* 


(2.41) 


Here  the  transpose  of  the  transpose  of  (AA")(AA")*  remains 
unchanged.  Hence  the  left  hand  side  of  Eq  (2.40)  is  symmetric  and 
consequently  so  is  the  right  hand  side  of  this  equation. 

Next,  assume  that  MT*=0  holds  by  hypothesis.  By  making  use 
of  theorem  2-2  which  guarantees  the  existence  of  a  solution  x  to 
the  equation  Ax=b  when  MT*=0  (MT*=0  was  the  consistency 
condition  for  a  general  solution  of  solution  ^  +  Nz  ). 

Apply  the  results  of  this  theorem  to  the  columns  of  T*.  Let  b  be  the 
columns  of  T*  where  b  may  vary  as  only  one  to  all  the  columns  of 
T*.  Then  by  theorem  2-2,  there  exists  a  matrix  X  such  that  AX=b 
holds  where  b  ranges  over  the  columns  of  T*,  A  is  fixed,  and  X 
arises  from  each  solution  of  AX=b  .  Now  T*=AX  implies  that 


((T)T  =  T  =  (AX)*  =  X*A‘ 


(2.42) 


Since  A~  =  ST 


A’  =  S(X*A*) 
then  for  Z  =  SX* 
(SX*)A*  =  ZA* 


Finally  it  follows  from  the  earlier  equation  that 

(A  A’lA  A‘)*  =  A  A'(A‘)*A* 

=  A(A'(A‘)*A* 

=  A(STA‘*A*  ) 

=  ASX*A*A‘*A* 

=  AZ(A*A'*A*) 

=  AZ(A  A'A)* 

=  AZA* 

=  A(ZA*) 

=  AA' 


So  from  Eqs  (2.42)  through  (2.45),  equality  holds  for 
(A  A'XA  A~)*KA  A‘) 

and  by  this  equality  symmetry  also  holds. 


Theorem  2i4.  Recall  the  definition  of  the  nonsigular  matrices  R 


and  C 


R  =  \  and  C  =  [S  M] 


and  the  result  of  theorem  2-1,  namely: 


(2.9) 


RaC  r 

C  0 


If  I  0  T 
0  0  IVI 
S  I N  0 


(2.10) 


Then  the  equation  Ax  =  b  has  a  solution  if  and  only  if  N*b  =  0  and 
the  general  solution  of  Ax  =  b  is  given  by 


^general  solution"^  S  b  +  M  z,  Vz 


(2.47) 


Proof.  Recalling  the  matrix  multiplication  which  lead  to  the 
form  of  Eq  (2.10)  in  theorem  2-1,  now,  the  transpose  of  this  sequence 
of  multiplications  is  taken  prior  to  the  actual  multiplication.  Thus 
from  the  original  multiplication  sequence 


R  U  A  I  0  _  RA  R  C  u  _  RAl  F: 

0  I  I  0  0  1  I  0  0  I  C  0 


and  by  taking  the  transpose  of  Eq  (2.8)  the  following  is  obtained 


* 


*  *•  b'* 


ifo]  [a  i]  [r*o]  _[c*a  c*]  [r*o]_[c*ar*c*' 
0  ij  [l  oj  [o  lj“  [  I  0  J  [o  iJi  R*  0  . 


Extending  the  results  of  theorem  2-1  to  the  form  of  Eq  (2.48) 
provides  the  following  equality 


-o*  .-.*1 

lAR  l 
[  R  u  J 


Ir-lO  ^ 

=  cT^  n' 

T^llVll  0 


(2,49) 


where  now 

R*  =  [T*  M*] 


Since  R  and  C  are  nonsingular  matrices  by  hypothesis,  the  following 
series  of  implications  can  be  made 

A*x=b  has  a  solution  x 


if  and  only  if  C*A*x  =  C*b  has  a  solution  x 

if  and  only  if  C*A*R*y  =  C*b  has  a  solution  y  where  x  =R*y 

The  last  implication  can  be  rewritten  using  Eq  (2.49).  Thus  the  last 
implication  becomes 
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where  y  is  now  defined  as 


with  X  consequently  rewritten  as 


X  =  [T*  M*] 


w 


Now  the  last  implication  as  rewritten  in  Eq  (2.52)  can  be  simplified 
as  follows 

Ir  0 
0  0 


V  = 


ir 


(2,53) 


It-  0 

w 

0  0 

z 

[S^b  N"b] 


(2.54) 


[W  0]  =  [S*b  N*b] 


(2.55) 


For  the  simplification  in  Eq  (2.53)  through  Eq  (2.55)  to  hold,  the 
following  equalities  must  be  true 

W  =  S*b  and  N*b  =  0 


X  =  T  W  +  M  z 


Finally  using  the  new  definition  of  W  provides  the  desired  result 


^general  solution"  T  S  b  +  M  z,  Vz 


(2.56) 


Corollary  2.  The  columns  of  M*  form  a  basis  for  the  null  space 
of  A*,  referred  to  by 

Proof.  The  proof  for  corollary  2  is  the  similar  to  that  of 
corollary  1,  but  follows  from  the  representation  of  the  general 
solution  of  X  given  in  theorem  2-6. 


From  theorem  2-6  the  last  two  remaining  generalized  inverses 
can  be  obtained.  First,  using  the  results  established  in  theorem  2-6, 
the  Aj  2  4  generalized  inverse  will  be  obtained  in  the  following 
theorem.  The  final  generalized  inverse  of  interest  will  be  obtained 
using  the  definition  of  an  Aj  2  3  4  generalized  inverse  established 
earlier. 


Theorem  2-7  (12).  If  N*S=0  then  the  matrix  ST  is  an  Aj  2  4 
generalized  inverse  of  A. 


Proof.  To  prove  ST  is  an  Aj  2^4  generalized  inverse,  it  must  be 
shown  that  ST  satisfies  the  following  equations 


(2.57) 

^1,2,4  ^^1,2,4  "  -^1,2,4 

(2.58) 

^^1,2,4^^*  "  ^^1,2,4^^ 

(2.59) 

By  the  results  of  theorems  2-3  and  2-4  it  follows  that  Aj  2  4 
satisfies  Eqs  (2.57)  and  (2.58).  To  show  Eq  (2.59)  holds  it  must  be 
shown  that  (Aj2  4A)*  is  symmetric.  This  proceeds  in  a  manner 

similar  to  the  proof  of  theorem  2-5.  First  let  A”=ST  and  then 
show  Eq  (2.60)  is  symmetric. 

(A'A)*(A‘A)  =  (A'A)  (2.60) 

Using  the  definition  of  symmetry  mentioned  in  theorem  2-5,  the 
transpose  of  the  left  hand  side  of  Eq  (2.60)  is  taken  as  follows 

((A'A)*(A'A))*  =  (A'A)*((A'A)*)*  =  (A"A)*(A'A)  (2.61) 

However,  this  symmetry  it  does  not  establish  equality.  To  show 
equality  let  N*S=0,  and  then  by  using  the  results  of  theorem  2-6, 
there  must  exist  an  X  such  that  A*X=S.  Using  this  result  and  the 
previous  definition  of  A~  leads  to  the  following 


(A'A)nA'A)  =  (A*A~^)(A'A) 


=  (A*A‘*)(A*XTA) 

=  (A*A‘*A*)(XTA) 

=  (A*)(XTA) 

=  (A*XT)(A) 

=  (A‘A) 

Therefore,  the  ST  matrix  also  satisfies  Eq  (2.59)  and  is  thus  an  Aj  2  4 
generalized  inverse. 

■ 

The  final  inverse,  the  Aj  2  3  4  generalized  inverse  is  obtained 

from  theorems  5  and  7  using  the  defintions  illustrated  in  figure  2 
earlier.  Therefore,  if  MT*=0  and  N*S=0,  the  ST  matrix  must  satisfy 
■  Eqs  (2.1)  through  (2.4).  Then  by  definition,  the  matrix  ST  is  an 

Aj  2  3  4  generalized  inverse. 


Ill  The  ST  Computation 


Chapter  2  reviewed  the  theory  behind  the  ST  method  for 
computing  the  generalized  inverse  of  a  matrix,  and  briefly  introduced 
this  technique  through  an  example.  This  chapter  will  examine  the 
computer  implementation  of  this  technique  in  greater  detail  and  also 
discuss  an  interesting  adaptation  of  the  ST  method  which  will  make 
it  more  immediately  useful  and  available. 

In  respect  to  the  ST  method,  the  algorithm  behind  this 
computational  technique  will  be  examined  from  several  perspectives, 
including  time  and  space  complexities  plus  numerical  accuracy 
considerations.  To  rounded  out  the  discussion  and  put  it  in  the  proper 
perspective,  the  ST  method  will  be  compared  to  other  techniques 
currently  used  to  computed  various  generalized  inverses.  The  intent 
of  this  comparison  is  to  try  to  determine  what  possible  advantages  or 
potential  disadvantages  the  ST  method  may  possess  in  relation  to 
these  other  techniques. 

While  determining  the  potential  advantages  and  dis¬ 
advantages,  the  discussion  on  comparisons  will  also  serve  to  introduce 
the  adaptaion  of  the  ST  method  to  a  well  known  algorithm.  This 
adaptation  is  based  on  a  singular  value  decomposition  algorithm.  The 
ST  method  will  be  used  to  short  cut  a  method  used  for  computing  an 
A,  2  3  4  inverse  based  on  the  results  of  the  singular  value 


decomposition  algorithm.  This  discussion  will  conclude  the  chapter 
with  a  detailed  explanation  and  an  example  of  this  modified 


computational  technique  for  finding  an  AJ234  generalized  inverse. 


Ih£  SI  Method:  A  Synthgtid  Approach 

The  ST  method  can  be  described  as  a  five  step  technique  for 
finding  four  generalized  inverses  of  a  given  matrix.  Of  these  five 
steps,  two  steps  form  the  computational  heart  of  the  technique. 

These  steps,  step  2  and  step  4,  utilize  a  modified  Gauss- Jordan  and 
the  modified  Gram-Schmidt  orthogonalization  algorithms  respectively. 
A  brief  overview  of  each  step  will  be  given.  Following  the  overview, 
the  computations  surrounding  step  2,  the  modified  Gauss- Jordan 
algorithm,  and  step  4,  the  modified  Gram-Schmidt  orthogonalization 
algorithm,  will  be  examined  in  detail.  Based  primarily  on  steps  2 
and  4,  the  time  and  space  complexity  of  the  ST  method  will  be 
discussed.  Next  some  of  the  subtleties  involved  in  implementing  the 
ST  method  will  be  addressed  from  a  very  general  view.  Included  in 
this  discussion  will  be  numerical  accuracy  considerations.  Finally,  this 
technique  will  be  compared  to  others  used  to  compute  various 
generalized  inverses. 

Algorithm  Overview.  The  example  in  chapter  1  used  a  new 
representation  defined  by  theorem  1-1,  namely  Eq  (2.10),  as  the 
starting  point  for  computing  four  generalized  inverses  To  obtain  this 
representation,  step  !  begins  by  augmenting  the  matrix  of  concern. 


Figure  7.  Dimensions  of  Augment  Matrix  and  Identities 

say  matrix  A,  with  an  identity  matrix  beside  it  and  an  identity 
matrix  below  it.  If  the  given  matrix  A  is  of  dimension  mxn,  then  the 
identity  matrix  augmented  beside  A  must  be  of  dimension  mxm.  The 
identity  matrix  augmented  below  A  must  have  dimensions  nxn.  Thus 
the  result  of  step  1  is  a  square  matrix  with  dimensions  (m+n)  by 
(n+m)  as  depicted  in  figure  7.  This  matrix  representation  was  shown, 
by  theorem  2-1,  to  be  similar  to  the  representation  which  will  now 
be  obtained  in  step  2.  Before  leaving  step  1,  note  the  zero  matrix  in 
the  lower  right-hand  corner  of  figure  7.  While  it  is  not  involved  in 
the  actual  computation  and  is  never  changed,  it  has  a  potential  use 
which  will  be  pointed  out  later. 

Now  to  obtained  the  next  representation  of  the  augmented  A 
matrix,  step  2  utilizes  a  slightly  modified  Guass-Jordan  algorithm  for 
reducing  A  to  an  identity  matrix  with  dimensions  equal  to  the  rank 
of  A.  Applying  this  algorithm,  the  representation  of  figure  7  becomes 
the  famiilar  representation  illustrated  in  figure  8.  Important  to  note 
from  figure  8  is  how  the  subma trices  S,  T,  M,  and  N  arise.  First, 


0 


I  '!m-r)xri-i 

0  I  0  M 

fix  f  fiAx  ri~  t  / 

S  N  0 

Figure  8.  Partitioning  and  Dimensions  of  Submatrices 

these  submatrices  are  the  consequence  of  reducing  of  A  to  1^.  where  r 

is  the  rank  of  A.  Second,  the  partitioning  of  what  was  originally  the 
two  identity  matrices  is  determined  by  r.  Thus  the  respective 
dimensions  of  these  submatrices  are  S^^^,  and 

fjnx(n-r)  illustrated  in  figure  8.  As  mentioned  earlier,  this  repre¬ 
sentation  forms  the  foundation  for  computing  the  four  generalized 
inverses  of  interest  in  this  work,  and  its  validity  is  based  on  the 
theory  reviewed  in  chapter  2. 

In  the  representation  of  figure  8,  the  computation  of  the  first 
generalized  inverse  of  interest  follows  directly  in  step  3.  Here  the  A|2 

inverse  is  computed  simply  from  the  product  of  the  S  and  T 
submatrices. 

Computing  the  next  two  inverses,  Aj2  3  and  A.^2  4’ 

the  use  of  the  modified  Gram-Schmidt  orthogonalization  algorithm. 
Again  using  the  representation  of  figure  8,  each  row  in  T  is 
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orthogonalized  to  every  row  in  M.  Now  taking  the  product  of  the 
original  S  submatrix  and  the  new  T  submatrix  gives  the  Aj  2  3 

generalized  inverse.  This  process  is  defined  as  step  4A.  Step  4B  again 
uses  the  same  orthogonalization  process  but  now  it  is  applied  to  each 
column  of  S.  Each  column  of  the  submatrix  S  is  orthogonalized  to 
every  column  of  N.  The  product  of  the  new  S  and  old  T  submatrices 
produces  the  Aj  2  4  generalized  inverse.  Since  both  inverses  are 

directly  based  on  the  results  of  step  3,  the  two  separate  computations 
were  labeled  as  steps  4A  and  4B. 

The  final  inverse  the  AJ234  generalized  inverse,  follows  from 

both  steps  4A  and  4B.  This  step  is  naturally  labeled  step  5,  and  is 
composed  simply  of  the  product  of  the  orthogonalized  S  and  T 
submatrices. 

The  relationship  of  these  five  steps  along  with  their  important 
characteristics  is  summarized  on  the  following  page  in  figure  9. 

Modified  Guass-Jordan  Algorithm.  The  Guass-Jordan 
algorithm  is  well  documented  (22:71,  27:417,  24,  1,  17),  and  the  details 
of  this  method  will  not  be  discussed.  Instead,  the  aspects  of  this 
method  with  particular  importance  to  the  ST  method  will  be 
addressed.  In  particular,  this  includes  the  addition  of  a  column 
operation  to  reduce  a  given  matrix,  and  a  change  to  what  might 
normally  be  used  as  a  pivoting  scheme  for  reducing  the  given  matrix. 
These  two  issues  are  discussed  below  as  they  pertain  to  the  modi¬ 
fications  of  the  Guass-Jordan  algorithm  adopted  for  the  ST  method. 


STEP  2:  REDUCE  A 


Reduce  A  to  an  identity  matrix  with 
dimenaion  equal  to  the  rank  of  A  to  form: 


Ir|  0 

T‘ 

0  1  0 

M 

SIN 

0 

STEP  3:  FORM  A  i  ,2  INVERSE 


Form  Ai  2  from  the  product  of  the  sub- 
matricea  S  and  T.  At  ^2  satiafiea  the 
following  conditiona. 

A  Ai  2  A  -  A 

A|,2  AAi,2  =Ai,2 


STEP  4A:  FORM  A  1,2,3  INVERSE 


Form  Ai  2,3  from  the  product  of  the  aub- 
matriceaSand  T  where  each  row  of  Tia 
orthogonal! zed  to  the  rowa  of  the  aub- 
matrix  M.  Ai  ,2,3  aatiafiea  the  following 
conditiona; 

A  Ai  ,2,3  A  =  A 
Ai,2,3  a  Ai,2,3  =Ai,2,3 
<A  Ai,2,3  )*  =  (A  Ai  ,2,3  ) 


STEP  4B:  FORM  A  1,2,4  INVERSE 


Form  Ai  ,2,4  from  the  product  of  the  aub- 
matricea  S  and  T  where  each  column  of  S  ia 
orthogonalized  to  the  columna  of  the  aub- 
matrix  N  A  1,2,4  aatiafiea  the  following 
conditiona: 

AAi,2,4  a  =  A 
Ai,2,4  a  At  ,2,4  *Ai,2,4 
(Ai,2,4  a)*  =  (Ai,2,4  A) 


STEP  5  FORM  A  1,2, 3,4  INVERSE 


Form  Ai  2,3,4  from  the  product  of  S 
in  atep  4B'and  T  in  atep  4A,  Ai  ,2,3,4 
aatifiea  the  following  conditiona: 

A  Ai  ,2,3,4  A  =  A 
Ai  ,2,3,4  A  A,  ,2,3,4*  A,  ,2,3,4 

(A  At  ,2,3,4  )*  =  (A  At  ,2,3,4  ) 
(At  ,2,3,4  A)*  *  (At  ,2,3,4  A) 


Figure  10.  Underdetermined  Case 


Typically  the  Guass- Jordan  algorithm  is  used  to  reduce  a 
square  matrix  to  row-echelon  form.  In  this  application,  matrices  are 
not  square  and  row-echelon  form  is  not  always  sufficient.  Thus  the 
Guass- Jordan  method  was  modified  to  reduce  a  non  square  matrix  to 
the  form  described  in  step  2  or  likewise,  derived  in  Eq  (2.10).  In  some 
cases,  the  row-echelon  form  produced  by  a  typical  implementation  of 
the  Gauss-Jordan  method  would  be  sufficient  for  non  square 
matrices,  but  given  the  case  of  a  matrix  with  dimension  mxn  where 
n  is  greater  than  m,  the  desired  form  would  not  be  reached.  This 
underdetermined  case  is  illustrated  in  figure  10.  In  figure  10  part  of 
the  original  matrix  was  reduced  to  an  mxm  identity  matrix,  but  this 
left  the  remaining  (n-m)  columns  non-zero.  For  the  ST  method  to 
work,  the  remaining  columns  in  the  underdetermined  case  must  be 
reduced.  This  problem  is  eliminated  by  was  using  both  row  and 
column  operations  to  reduce  a  given  matrix  Row  operations  are  used 
to  first  reduce  the  matrix  to  an  upper  triangular  form  with  rank 
equal  to  that  of  the  original  matrix.  Next  applying  column  operations. 


the  remaining  elements  above  the  diagonal  are  eliminated.  This 
provides  the  desired  form.  While  it  would  have  been  possible  to  flag 
the  underdetermined  cases  for  special  treatment,  this  was  not  done. 
Instead,  in  the  interest  of  generality  and  numerical  accuracy 
considerations,  which  will  be  elaborated  on  later,  a  method  using 
both  row  and  column  operations  was  adopted. 

The  addition  of  column  operations  to  the  existing  Guass-Jordan 
algorithm  required  an  accopanying  change  to  the  pivoting  scheme 
normally  employed  with  this  algorithm.  Typically,  a  satisfactory 
pivoting  scheme  for  most  cases  would  simply  search  each  column 
below  the  diagonal  for  the  best  pivot  element  before  attempting  to 
eliminate  the  remaining  elements  (22:187-188;  24:38).  This  pivoting 
scheme  works  well  when  only  row  operations  are  used  to  reduce  a 
matrix,  but  ignores  the  subsequent  column  operations  which  are 
required  for  the  ST  method.  Since  column  operations  possess  the 
same  potential  problems  found  in  row  operations,  namely  subtracting 
numbers  of  similar  size,  any  pivoting  scheme  used  must  also  be  able 
to  look  ahead  for  the  subsequent  column  operations  (27:593).  One 
immediate  technique  would  be  adopting  a  complete  pivoting  scheme 
where  the  entire  matrix  is  searched  for  the  best  pivot  element.  The 
obvious  problem  with  this  approach  is  the  computational  expense  of 
searching  the  entire  matrix.  Since  complete  pivoting  schemes 
normally  only  provide  two  to  four  times  better  results  at  a  high  cost 
(24:39)  this  method  was  not  used.  Instead  a  compromise  was  designed 
which  uses  at  the  most  two  column  searches  and  a  single  row 
search.  Thus,  from  a  given  diagonal  element,  a  column  search  is 
conducted  below  the  diagonal  element  followed  by  a  row  search  to 
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the  right  of  the  diagonal  element.  If  the  pivot  element  was  chosen 
from  the  row  search,  a  third  search  is  conducted  on  the  new 
column.  This  scheme  is  employed  during  the  row  reduction  operations 
as  the  given  matrix  is  reduced  to  an  upper  triangular  form.  The 
intent  is  to  provide  a  good  computing  framework  for  the  following 
column  operations  when  the  form  of  the  matrix  is  less  flexible.  Here 
less  flexible  means  it  is  difficult  to  pivot  during  column  operations 
without  loosing  a  zero  element.  The  alternative  of  pivoting  during 
column  operations  with  loosing  zeroes  would  begin  to  approach  the 
computational  burden  of  complete  pivoting  since  the  matrix  would 
have  to  be  rearranged  back  to  the  form  of  Ij,  at  some  point.  Thus  this 

scheme  avoids  the  expense  of  complete  pivoting  but  at  the  same  time 
provides  a  good  basis  for  the  column  operations. 

Spreading  the  numerical  operations  between  both  the  rows 
and  columns  has  the  added  advantage  of  reducing  the  number  of 
operations  performed  on  the  elements  in  each  of  the  two  augmented 
identity  matrices.  Instead  of  concentrating  the  changes  on  a  single 
identity  matrix  like  the  unmodified  Guass-Jordan  method  does,  this 
reduction  technique  spreads  the  arithmetic  operations  across  two 
identity  matrices  in  an  almost  equal  fashion.  The  result  is  less  opera¬ 
tions  on  individual  elements  of  the  submatrices  S,  T,  M,  and  N.  Since 
repeated  simple  arithmetic  operations  accumulate  round-off  error, 
and  subsequent  round-off  errors  are  the  result  of  the  current  error, 
the  error  in  each  submatrice  should  accumulate  more  slowly(27:593). 

In  summary,  one  change,  column  operations,  was  added  to 
the  Guass-Jordan  algorithm  to  adapt  it  for  the  ST  method.  An 
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accompanying  pivoting  scheme  was  developed  to  account  for  the 
column  operations.  The  combined  result  is  an  effective  routine  which 
is,  at  the  least,  as  numerically  accurate  as  the  well  published  version 
of  the  Guass- Jordan  method. 


Modified  Gram-Schmidt  Orthoeonalization  Algorithm.  Steps  4A 
and  4B  require,  depending  on  the  rank  of  the  matrix  of  interest, 
specific  rows  of  the  T  and  M  matrices  and  columns  of  the  S  and  N 
matrices  to  be  orthogonal.  To  orthogonalize  the  required  rows  or 
columns,  the  modified  Gram-Schimdt  orthogonalization  process  was 
used.  This  is  a  modified  version  of  the  original  Gram-Schmidt 
algorithm  and  differs  only  in  the  order  calculations  are 
performed  (24:151). 

The  use  of  the  modified  Gram-Schmidt,  in  the  ST  method, 
differs  a  little  from  its  typical  description  since  the  orthogonalization 
process  depends  on  the  rank  of  the  matrix  of  interest.  Thus  algorithm 
is  first  described  below  in  its  most  general  representation.  Its 
application  to  the  ST  method  is  then  explained  and  demonstrated 
with  an  example. 

The  modified  Gram-Schmidt  orthogonalization  algorithm  is 
described  with  the  following  short  segment  of  a  pseudocode.  In  this 
desciption,  the  purpose  is  to  orthogonalize  each  vector  denoted  with 
the  small  letter  a  and  subscripts.  The  algorithm  is  as  follows: 
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Given  the  vectors  aj  for 

For  k=l  to  M  do 
For  j=k+l  to  M  do 
aj  -  aj-[(aj.a^)/(a^*a^)Ja^ 

End  of  j  loop 

End  of  k  loop 

The  algorithm  description  also  serves  to  point  out  two  useful  features 
of  this  algorithm,  its  use  of  storage  and  its  simple  and  direct  nature. 

In  respect  to  storage,  the  modified  vectors  replace  themselves  so  no 
temporary  storage  is  needed  for  intermediate  results  which  reduces 
its  storage  over  the  original  Gram-Schmidt  process.  Second,  the 
process  is  completely  described  in  a  few  short  lines  which  attests  to 
its  simple  and  direct  nature. 

Since  the  ST  method  requires  that  each  individual  row  of  the 
T  submatrix  be  orthogonal  to  the  M  submatrix,  the  change  to  the 
basic  algorithm  described  above  involves  only  the  loop  indices.  For  the 
ST  method  and  the  T  and  M  submatrices,  the  outer  loop  is  limited  by 
the  rank  of  the  original  matrix.  The  inner  loop  is  controlled  by  the 
number  of  rows  in  the  augment  matrix  A.  In  terms  of  the  algorithm, 
this  restricts  the  actual  number  of  arithmetic  operations  performed 
over  the  case  where  a  set  of  vectors  or  rows  for  matrices  are  fully 
orthogonalized  For  the  case  of  the  T  and  M  submatrices,  each  row  in 
T  is  fully  orthogonalized  to  the  rows  of  M,  but  the  rows  of  T  are  not 


orthogonalized  to  each  other.  The  same  applies  to  the  orthogonaliza- 
tion  of  each  column  of  S  to  the  entire  N  submatrix.  The  following 
example  demonstrates  the  orthogonalization  process  used  for  the  ST 
method. 


Example  2.  Given  the  following  reduced  matrix  with 
dimension  6x4  and  a  rank  of  2; 
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Orthogonalize  each  individual  row  in  T  to  all  the  rows  of  M.  Here  the 
orthogonalization  process  begins  with  the  bottom  row  in  M  so  the 
identity  matrix,  Ij,  will  be  preserved.  If  the  process  were  started  with 
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the  first  row  in  T  and  worked  down,  the  Ij.  form  obtained  in  the 

previous  step  would  be  lost.  This  is  because  oper^.  ions  are  not  per¬ 
formed  on  only  the  rows  of  T  or  M,  but  on  the  entire  matrix. 
Therefore  care  must  be  taken  to  preserve  the  work  done  to  this 
point.  Thus  with  the  outer  loop  index  set  to  6  (the  fourth  row  in  M), 
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and  the  inner  loop  is  initialized  to  5,  the  algorithm  iterates  down  to  1 
This  results  in  the  following  changes  to  the  T  and  M  submatrices 
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After  the  first  iteration  of  the  outer  loop,  the  last  row  in  M  is 
orthogonal  to  each  row  above  it  including  the  two  rows  in  T.  The 
next  pass  of  the  outer  loop  will  begin  at  the  third  row  of  M  and 
orthogonalizes  it  to  each  row  above  it.  The  resulting  matrix  is: 
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The  process  is  repeated  again  for  the  second  row  of  M  orthogonalizing 
it  to  each  row  above  it  as  the  inner  loop  iterates  from  3  to  1.  The 


resulting  matrix  after  this  iteration  is 
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0  0  4  “  0 . 0  o  0 . 1 Z  -0 . 1 2  0 . 0 0  0 . 0  4  I 


0,67  1.25-0.37  0.67  0  00-0.67 

0  25  0.50-0  25  0.25  100-0.25 
0  20-0.40  0.60  1  00  0  00-0  20 
0.50  1.00  1.00  0  00  0  00  0  50 

1  0  0  0  0 0  0  0 0  0  0 0  0  0 0  1  0 0 


Finally,  the  first  row  in  M  is  orthogonalized  to  each  row  above  it, 
namely  the  two  rows  in  T.  The  orthogonalizauon  process  is  stopped  at 
this  point  by  the  rank  of  the  original  matrix.  The  final  T  and  M 
submatrices  are 
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The  same  process  is  performed  for  the  submatrices  S  and  N  but  on  a 
column  basis.  This  time  the  algorithm  only  completes  two  outer 
iterations  to  meet  the  orthogonalization  requirement.  Again  the 
orthogonalization  process  is  stopped  by  the  rank  of  the  original 
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matrix.  The  final  matrix  with  the  results  of  both  orthogonalization 
processes  and  the  unchanged  identity  matrix  is  shown  below.  Step  5 
follows  directly  from  this  final  form  as  the  product  of  the  orthogon- 
alized  S  and  T  submatrices 
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One  remaining  aspect  of  this  technique  not  yet  discussed  is  it’s 
numerical  accuracy.  As  described  above,  the  modified  Gram-Schmidt 
process  is  a  numerically  stable  technique  (24:152).  This  is  based  on 
completely  orthogonalizing  every  vector  or  row  in  this  case  to  every 
other  vector  or  row.  In  the  ST  method,  the  orthogonalization  process 
is  not  complete,  but  instead  limited  by  rank.  In  terms  of  the  arith¬ 
metic  operations  performed,  this  means  the  modified  Gram-Schmidt, 
as  adapted  for  the  ST  methods,  performs  less  arithmetic  operations 
than  the  general  case  which  is  described  above  as  numerically  stable. 
Thus  the  accuracy  of  the  Gram-Schmidt  orthogonalization  possess  in 
this  application  should  be,  at  the  least,  as  good  as  the  above  reference 
suggests. 


Time  and  Space  Complexities.  Two  measures  of  any 
algorithm  are  its  speed  and  its  storage  requirements,  and  these 
measures  are  often  referred  to  as  an  algorithm's  time  and  space 
complexities  respectively  (25:24-39).  For  the  ST  method,  the  speed  or 
time  complexity  analysis  will  be  limited  to  the  two  dominating 
processes,  the  modified  Guass-Jordan  algorithm  and  the  modified 
Gram-Schmidt  orthogonalization  algorithm.  For  storage  requirements 
or  the  space  com.plexity,  only  the  augmented  matrix  will  be 
considered.  These  limitations  are  imposed  to  make  the  analysis  more 
general  and  consequently  more  accurately  compared  to  other  existing 
algorithms. 

The  time  complexity  of  the  ST  method  is  based  on  counting 
numerical  operations  performed  in  the  two  dominating  processes.  This 
is  even  further  limited  by  using  only  floating-point  operations  as  is 
commonly  done  (17:78). 

Thus  the  basis  of  the  modified  Guass-Jordan  reduction 
becomes  two,  three  level  loops.  Summarized  by  the  row  operations 
performed  on  a  matrix  with  dimensions  mxn,  the  loops  consist  of 

For  k=l  to  M 
For  j=k+l  to  M 
For  i=k+l  to  (M+N) 


For  column  operations  the  loops  are  identical  with  the  exception  of 
the  outer  loop  limits.  The  limit  for  the  outer  two  loops  is  N.  Using  a 
common  approximation  for  the  above,  the  total  row  and  column 
loops  are  equated  to  m^/3  or  n^/3  floating-point  operations 
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respectively  (17:80).  The  total  operation  count  for  the  matrix 
reduction  is  then  simply  the  sum  of  each  which  is  (m^+n^)/3. 

For  the  modified  Gram-Schimdt  orthogonalization  process,  the 
time  analysis  is  the  same.  Each  orthogonalization  process  contains  a 
three-deep  nested  loop  as  illustrated  above  for  the  row  reduction 
operations.  In  the  worst  case,  with  a  matrix  rank  of  one,  the 
operation  count  for  the  row  orthogonalization  is  n^/3,  and  for  the 
column  orthogonalization  the  count  is  m^/3.  Again  the  combined 
count  is  simply  the  sum  of  each  which  is  (m^+n^)/3. 

Combining  the  operation  count  for  the  Guass-Jordan  and 
Gram-Schmidt  provides  an  overall  operation  count  of  (2/3)(m^+n^). 
This  applies  to  the  computation  of  an  Aj  2  3  4  generalized  inverse. 
Computing  a  lower  inverse  would  require  less  operations  down  to  a 
minimum  of  (m^+n^)/3  for  an  Aj  2  generalized  inverse.  No  other 

operations  performed  by  the  ST  method  would  appreciably  change 
the  cubic  derived  above.  While  including  these  other  operations  might 
change  the  coefficients  or  add  a  lower  term  to  the  overall  operation 
count,  the  cubic  nature  of  the  two  dominating  processes  will  still 
prevail  as  the  methods  chief  computational  burden  and  also  serve  as 
a  good  comparative  reference. 

The  space  complexity  is  based  on  the  number  of  memory 
locations  occupied  by  the  augment  matrix.  While  the  actual  imple¬ 
mentation  used  considerably  more  (extra  storage  was  used  for 
validation  and  demonstration  purposes),  the  heart  of  the  algorithm 
requires  only  enough  storage  for  the  original  matrix,  two  identity 
matrices,  and  the  resulting  generalized  inverse.  The  analysis  flows 


very  directly  from  this  perspective.  For  the  original  matrix,  a  block 
of  memory  determined  by  the  dimensions  of  the  matrix  is  needed.  If 
given  a  matrix  with  dimensions  m  by  n,  then  m  times  n  memory 
locations  are  needed.  Since  the  original  matrix  is  augmented  by  two 
square  identity  matrices,  an  additional  m^+n^  locations  are  needed. 
The  resulting  generalized  inverse  needs  n  times  m  locations  since  it 
has  dimensions  n  by  m  based  on  the  original  matrix.  Combining  the 


space  requirement  for  the  original  matrix,  the  two  identity  matrices, 
and  the  generalized  inverses  results  in  (m+n)^  locations. 

Other  Algorithm  md  Computation  Considerations.  Lumped 
under  this  section  are  the  subtleties  and  hidden  insights  enclosed  in 
the  computer  implementation  of  the  ST  method  which  might  not  be 
obvious  from  the  desciption  above  or  its  documentation.  Each  is 
addressed  individually  and  with  no  specific  order  in  the  paragraphs 
which  follow.  These  should  be  helpful  to  any  future  implementations 
of  this  method  for  computing  the  generalized  inverse  of  a  matrix. 

The  simplest  data  structure  for  implementing  the  augmented 
matrix  is  a  square  two  dimension  array.  For  a  matrix  with 
dimensions  mxn,  such  an  array  would  be  dimensioned  (m+n)  by 
(m+n).  Although  this  representation  is  the  simplest,  it  also  poses  one 
problem  regarding  storage.  In  the  space  analysis  the  complexity  was 
determined  to  be  (m+n)^.  Using  this  scheme  the  original  matrix  and 

the  two  identities  matrices  will  require  (m+n)^  locations  alone.  To 
keep  the  storage  to  the  minimum  found  in  the  space  analysis,  the 
computed  generalized  inverse  could  be  stored  in  the  unused  lower 
right  corner  of  the  (m+n)  by  (m+n)  array.  This  was  eluded  to  earlier 


(on  page  A4)  when  discussing  the  composition  of  the  augmented 
matrix  after  it  was  reduced.  While  this  may  require  a  little  extra 
effort  during  the  final  phase  of  the  program  to  place  the  product  of  S 
and  T  back  into  the  lower  right  corner  of  the  large  array,  the  trouble 
may  be  offset  by  the  ease  of  manipulating  only  one  data  structure 
and  still  using  minimum  storage.  A  second  alternative  would  be  to 
store  the  computed  generalized  inverse  on  top  of  the  identity  matrix 
which  was  denoted  as  Ij,.  There  are  also  problems  with  this  approach, 

namely  taking  the  transpose  of  the  inverse  while  storing  and  retriev¬ 
ing  it  from  this  area,  but  again  the  convenience  may  out  weight  the 
trouble  of  manipulating  only  one  data  structure. 

The  use  of  a  conditional  statement  with  a  floating-point 
representation  introduces  a  problem  when  the  condition  is  testing  for 
equality.  In  the  ST  method,  this  problem  arose  during  the  matrix 
reduction  process  of  step  2  when  testing  for  zeroes.  Since  a  given 
matrix  element  may  never  exactly  equal  a  floating-point  zero,  due  to 
small  numerical  errors,  the  question  must  be  asked  how  should  the  a 
condition  be  setup  to  test  for  a  floating-point  zero?  One  method 
represents  a  zero  relative  to  the  scale  of  the  machine  used,  and  then 
bases  the  condition  on  the  absolute  value  of  this  relative  zero  (24:43). 
Taking  the  machine  used  for  the  original  development  which  had  16 
decimal  digits  of  accuracy,  a  zero  was  represented  as  5x10“^^.  Thus 
the  absolute  value  of  any  matrix  element  less  than  5x10"^^  was 
considered  zero.  Use  of  a  different  machine  would  require  an 
adjustment  to  this  representation  of  zero.  Additionally,  different  size 
matrices  would  require  adjustments  to  the  amount  of  error  allowed 
in  a  zero 
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There  are  certainly  few  if  any  applications  where  all  four 
generalized  inverses,  discussed  in  this  work,  are  needed  together. 

Thus  only  certain  parts  of  the  ST  method,  as  implemented  in 
appendix  B,  are  needed  for  a  given  application.  This  reduces  the 
complexity  of  a  potential  program  significantly  from  what  is  included 
in  appendix  B.  As  mentioned  above  under  the  time  complexity 
analysis,  this  results  in  a  range  of  operation  counts  from  (l/3)(m^+n^) 
to  (2/3)m^+n^)  floating  point  operations.  The  same  applies  to  the  other 
routines  included  in  appendix  B  supporting  the  computation  of  a 
given  generalized  inverse.  These  include  the  multiplication,  printing, 
and  input  routines.  In  any  given  application,  these  routines  could  be 
reduced  significantly  to  support  the  given  application's  requirements 
lessening  the  computational  burden  of  the  ST  method. 

Another  aspect  of  applying  the  ST  method  concerns  the 
environment  where  the  application  would  actually  reside.  Consider 
an  environment  with  a  large  mathematical  software  package.  Such  a 
package  could  be  used  to  piece  together  an  effective  ST  implement¬ 
ation  with  minimal  modifications.  For  example,  by  calling  a 
Gauss-Jordan  algorithm  followed  by  any  necessary  column  reductions 
and  a  call  to  a  matrix  multiplication  routine,  an  Aj  2  inverse  could  be 

computed  with  little  effort.  There  are  certainly  many  other  possibili¬ 
ties  which  are  made  possible  due  to  the  modular  nature  of  the  ST 
method. 

One  final  topic  concerns  algorithm  validation.  Two  methods 
were  used  to  validate  this  effort.  The  most  useful  was  testing  the 
computed  inverse  against  the  one  to  four  conditions  it  is  suppose  to 
satisfy.  These  conditions,  Eqs  (2-1)  through  (2-4),  provide  a  reliable 
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means  of  testing  the  algorithm.  The  second  method  was  simply  to 
compare  the  computed  results  to  a  known  correct  result.  Appendix  A 
provides  a  detailed  example  showing  the  use  of  both  of  these 
methods.  The  conditions  are  tested  as  each  inverse  is  computed,  and 
the  final  unique  Aj  2  3  4  generalized  inverse  matches  the  known 

solution  found  in  the  referenced  text.  Of  the  two  methods  used,  the 
most  useful  and  informative  was  using  the  four  conditions.  This 
method  also  provides  insight  into  possible  numerical  inaccuracies 
since  the  results  can  be  compared  to  the  correct  solution.  This  is  the 
same  as  taking  the  classical  inverse  and  then  multiplying  it  and  the 
original  inverse  together.  The  resulting  identity  gives  an  indication  of 
how  accurate  the  computed  inverse  is  just  as  using  the  four  Penrose 
conditions  does  for  a  generalized  inverse. 

Algorithm  Comparison.  When  comparing  two  or  more 
(Algorithms  which  perform  the  same  basic  function,  but  do  so  by 
different  methods,  it  can  be  difficult  to  find  a  common  basis  for 
comparison  This  especially  true  in  this  case  where  there  are  many 
methods  for  computing  a  generalized  inverse.  Here,  the  comparison 
process  is  particularly  complicated  by  the  wide  variety  of  notation 
used  to  describe  such  computational  techniques.  Additionally  each 
description  is  normally  only  developed  to  provide  an  overview  of  the 
method  versus  the  detail  needed  for  an  indepth  comparison.  To 
overcome  these  difficulties,  six  general  categories  were  chosen  to 
compare  algorithms  Each  category  is  described  below 

The  first  category  deals  with  the  initial  representation  of  the 
matrix  This  category  was  chosen  since  the  process  of  augmenting  the 


initial  matrix  is  central  to  the  ST  method.  Thus  any  comparison 
involving  the  ST  method  should  include  this  particular  aspect.  From  a 
survey  of  several  different  methods,  the  second  category  was  chosen 
as  factorization  methods.  This  category  was  chosen  since,  in  one  way 
or  another,  all  algorithms  examined  performed  some  factorization  of 


!  Method 

1 

o 

3 

1 

i 

Characteristic 

ST 

(21 :142) 

...  ^+.-.t 

1 1.  y  II. 

(27:135) 

Intiai  Ivlatrir-: 

Fi  e  p  r  0  s  e  n  t  a  t  i  o  n 

El 

[A^'A  1  I] 

[A] 

Factorisation 

ROW  AMD 

COLUMN  OPS. 

ROW 

OPERATIONS 

SINGULAR  VALUE 

DECOMPOSITION 

Logical  F!ov/ 

DIRECT 

DIRECT 

MULTIPLE 

DECISIONS 

jHQSjyillllll 

4.'''.3(n)^ 

>  2.'''3(m''+n'‘) 

Space  CornpleKitsf 

2ri^ 

(rn+n)'^ 

j  Flexibility 

IHIH 

Restricted  by- 
rank 

Computes 

1  inverse 
_ 

Figure  11.  Algorithm  Comparison 


the  initial  matrix  in  the  process  of  computing  a  generalized  inverse. 
The  third  category  deals  with  an  algorithm's  general  structure. 
Here  general  structure  means  the  logical  flow  of  the  algorithm.  In 
other  words,  does  the  algorithm  require  extensive  logic  to  piece 
together  the  individual  parts  of  the  computation  or  does  it  flow 
directly  from  start  to  finish.  The  fourth  and  fifth  categories  are 


respectively  the  time  and  space  complexities  of  an  algorithm.  The 
final  category  deals  with  the  flexibility  of  the  algorithm.  In  this 
regard  flexibilities  concern  what  range  of  problems  can  the  method 
handle  or  what  constraints  are  levied  on  the  problem's  input. 

Using  these  six  categories,  the  ST  method  is  compared  to  two 
representative  methods  for  computing  a  generalized  inverse.  The 
comparison  is  summarized  in  figure  11  (Each  method  is  defined  and 
referenced  to  where  it  was  obtained  in  the  figure).  The  figure  shows 
a  number  of  general  results.  First,  the  ST  method  has  the  most 
complex  initial  representation,  but  this  representation  appears  to  lead 
to  an  advantage  in  the  other  characteristics  especially  in  flexibility. 
The  ST  method  is  the  only  one  of  the  three  which  can  compute  four 
generalized  inverses  without  any  restrictions.  Method  2  seems  to  be  a 
best  in  terms  of  storage,  computations,  and  factorization,  but  is 
restricted  by  the  rank  and  dimensions  of  a  given  matrix.  Method  3 
begins  with  the  simplest  initial  representation  but  quickly  becomes 
very  complex.  This  algorithm  will  be  examined  in  greater  detail  in 
the  discussion  of  the  adaptation  for  the  ST  method. 

In  summary,  it  hard  to  say  one  method  is  absolutely  better 
than  the  other,  but  certainly  the  ST  method  is  by  far  the  most 
flexible.  It  can  be  used  in  more  situations  than  the  other  two  and 
does  not  suffer  significantly  in  terms  of  the  other  categories  of 
comparison  due  to  this  advantage  in  flexibility.  Ignoring  all  other 
possibilities  and  using  the  six  categories  above,  the  ST  method  would 
probably  offer  the  best  choice  for  computing  a  given  generalized 


inverse. 


Ihs  21  Method  Adapted  lo  Singular  Value 


The  ST  method  will  be  adapted  to  the  results  of  a  singular 
value  decomposition  algorithm  to  obtain  an  Aj  2  3  4  generalized 

inverse.  This  will  begin  by  first  reviewing  an  algorithm  for  singular 
value  decomposition  to  provide  the  proper  perspective  for  the  ST 
adaption.  Next  the  representation  developed  for  the  ST  method  will 
be  derived  from  the  results  of  the  singular  value  decomposition 
algorithm  to  compute  the  Aj  2^3  4  generalized  inverse.  Finally  this 

approach  will  be  compared  to  another  method  which  uses  the  results 
of  a  singular  value  decomposition  algorithm  to  directly  compute  an 
Ai  2  3  4  generalized  inverse. 

Ih£  Singular  Value  Decomposition  Algorithm.  Given  a  matrix 
^mxn,  singular  value  decomposition  factors  the  matrix  into  the 
product  of  two  orthonormal  matrices  and  a  diagonal  matrix.  This 
product  is  defined  as 

,  -  ^  .-  .t 

A  =  L  Ur. 

To  understand  how  the  ST  method  is  adapted  to  the  results  of  the 
singular  value  decomposition,  the  Qj,  Q2,  and  2  matrices,  as  shown 

above,  will  be  developed  with  a  simple  example.  Then  the  S  and  T 
submatrices  will  be  derived  from  these  three  matrices.  The  process 


begins  by  first  choosing  the  matrix  A  as 


Next,  find  the  Q2  matrix.  This  is  begun  by  first  forming  A^A  and  then 

computing  I  A^A-M  I  with  the  goal  of  finding  the  characteristic  root 
of  A^A  as  follows 

=  [34] 

\.4  A-MI  =  !25-K|  =  0 
=  25 

With  a  characteristic  root  of  K  =  25,  the  process  continues  by  finding 
the  associated  characteristic  vector 


Using  the  characteristic  vector,  the  Q2  matrix  is  defined  as  follows 

-:i,“  !ij 

The  next  step  is  to  form  Qj.  This  begins  by  forming  AA^  and  then 


computing  I  AA^-M  I  to  finds  its  characteristic  roots  similar  to  what 
was  done  to  find  Q2.  In  this  case  the  process  produces  the  following 
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Notice  at  this  point  that  \  =  25  is  a  common  solution  for  both 
I  AA^-M  I  =  0  and  I  A^A-M  1  =  0.  The  next  step  uses  this  common 
solution  to  obtain  the  characteristic  vector  for  AA^  as  shown  below 


which  is  solved  for  x  as  the  following  pair  of  equations 

9xj  +  12x2  =  25xj 
12x,  +  16xo  =  25xo 


which  becomes 


-16xj  +  12x2  =  0 
12xj  -  9x2  ^  ^ 


Reduced,  these  equations  simply  become  Xj  =  (3/4)x2.  Then  by 
choosing  X2  as  4,  Xj  is  equal  to  3  and  the  characteristic  vector 
becomes 

r  ^  7 

4 

Normalizing  x  provides  the  following 


i  J 

Now  Q1  is  constructed  using  the  normalized  x.  The  first  column  of  Q1 
IS  X  leaving  one  additional  column  to  be  filled  to  make  the  product  Q1 
and  1  conform  This  second  column  is  filled  with  a  vector  making  Q1 
an  orthogonal  matrix.  In  this  case,  QI  is  defined  below  showing  this 
choice 
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The  last  matrix  to  be  found  by  singular  value  decomposition  is  2.  This 
matrix  is  the  product  of  Ql,  A  and  Q2  as  shown  below 


?  4 

I  -r-  F  f-,1  r^i 

-  I  ^  ;  [1] 

-  4  3  4  0 


The  diagonal  elements  of  1  display  the  singular  values  for  the  original 
matrix  A.  In  this  case,  the  decomposition  process  provides  only  one 
singular  value  of  A  which  is  simply  5.  This  will  become  important  as 
the  S  and  T  matrices  are  obtained. 

Deriving  ths.  S  Mid  I  Submatrices.  With  the  Ql,  1,  and  Q2 
matrices,  the  S  and  T  matrices  of  the  ST  method  can  now  be 
obtained  in  a  very  direct  manner.  First  a  matrix  is  formed  using  the 
Ql,  1  and  Q2  matrices  and  partitioned  as  shown  below. 
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This  representation  should  look  very  similar  to  the  representation 
developed  in  the  ST  method.  In  fact,  by  only  dividing  each  diagonal 
element  in  the  2  matrix  using  row  operations,  the  ST  method  is 
obtained.  This  is  equivalent  to  saying  divide  each  row  of  Q1  by  the 
corresponding  singular  values.  For  this  example,  there  is  only  one 
diagonal  element  or  singular  value  in  2  so  the  resulting  matrix  in  the 

ST  representation  becomes 
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As  described  many  times  before,  the  Aj  2,3,4  generalized  inverse 
simply  becomes  the  product  of  S  and  T.  For  this  example,  Aj  2,3,4 


This  adaptation  of  the  ST  method  is  certainly  simpler  than  the  more 
publisized  technique  defined  by 


..+  -t 


where 


To  avoid  the  pitfalls  associated  with  the  computation  of  a 
generalized  inverse,  potential  applications  were  modified  to 
circumvent  matrices  not  suitable  to  the  classical  method  of  matrix 
inversion.  Further,  some  potential  applications  probably  never  even 
considered  the  use  of  a  generalized  inverse,  especially  the  less  known 
generalized  inverses  like  the  Aj  2  3  •  One  such  application,  which 
probably  falls  somewhere  between  these  two  possibilities,  will  be 
adapted  to  use  an  Aj  2  generalized  inverse.  This  application,  the 

Lyapunov  matrix  equation,  will  be  solved  using  an  A|  2  generalized 
inverse.  The  solution  will  be  developed  in  detail  in  order  to  show  how 
the  Aj  2  generalized  inverse  fits  this  application. 


Jhs.  Lvapunov  Matrix  Equation 

The  necessary  and  sufficient  conditions  needed  to  show  the 
Lyapunov  equation,  Eq  (4.1),  has  a  solution  will  be  established.  This 

AX-XB>C  =  0  (4.1) 


where 


A  .  /-mxm 


solution  will  be  defined  in  terms  of  an  Aj  ^  generalized  inverse.  Two 

steps  will  be  required  to  accomplish  this  goal.  First,  the  solution  of 
the  Lyapunov  matrix  equation  must  be  defined  in  terms  of  two  pairs 
of  equations.  Each  pair  of  equations  will  trap  the  solution  of  the 
Lyapunov  matrix  equation,  if  a  solution  exists.  Establishing  each  pair 
of  equations  will  require  the  use  of  several  identities  which  will  be 
developed  in  the  bulk  of  this  section.  Second,  once  the  two  pair  of 
equations  are  established,  their  solution,  and  consequently  the 
solution  of  the  Lyapunov  equation,  will  be  expressed  using  an  Aj  2 

generalized  inverse.  Thus  a  general  common  solution  of  the  pairs  of 
equations  will  be  defined  to  express  the  solution  of  the  Lyapunov 
matrix  equation. 

Theorem  4-1.  Let  Eq  (4-1)  have  a  solution  X,  then  the  following 
matrices  are  simialr 
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Proof.  Let  X  be  a  solution  of  Eq  (4.1),  then  the  following  holds 
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Using  the  given,  namely  that  Eq  (4.1)  has  a  solution,  the  two  matrices 
of  Eq  (4.2)  are  similar.  Note  this  holds  for  any  X  since  the  matrices 
which  appear  in  the  left  hand  side  of  Eq  (4.3)  as  shown  below  in  Eq 
(4.6)  are  inverses 
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With  theorem  4-1,  it  has  been  shown  that  there  exists  an  X 
such  that  the  matrices  of  Eq  (4.2)  are  similar  and  consequently  Eq 
(4,1)  has  a  solution.  The  next  step  is  to  extend  this  result.  To  do  this, 
define  the  matrix  R  as  follows 


then 
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Next  find  the  determinant  of  [R-M]  as  follows 


(4.10) 


This  shows  that  I  R-M  I  factors  into  a  pair  of  polynomials,  namely 
the  characteristic  polynomials  of  the  matrices  A  and  B.  The  degree  of 
is  less  than  or  equal  to  m  and  the  degree  of  gg(\)  is  less  than  or 

equal  to  n.  Also,  fp,(R)  is  equal  to  zero  since  every  matrix  satisfies  its 

own  characteristic  equation. 

Now  using  the  results  of  theorem  4-1,  namely 
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let  tpp\)  be  defined  as  follows  for  ^  replaced  by  R 


where  U,  V,  M,  and  N  are  polynomials  in  the  matrices  A,  B,  and  C 
Any  polynomial,  f^(R),  commutes  with  R  where  the  coefficients  are 
complex  or  real,  this  results  in  the  following 
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This  implies  the  following  identities 


AU+CV  =  UA 
BV  =  VA 

AM+CN  =  UC+MB 
BN  =  VC+NB 


(4,15) 
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Again  recalling  the  results  of  theorem  4-1  in  the  form  of  Eq  (4.11) 
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This  equation  holds  for  any  solution  X  of  Eq  (4.1) 

AX-XB+C  =  0 

Thus  for  any  polynomial  f^(^)  given  by  Eq  (4.12) 
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and  under  the  similarity  transformation  as  defined  and  proved 
earlier  for  Eq  (4.2),  it  follows 
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where  f^(B)  is  not  equal  to  zero  and  *  is  not  of  concern.  Next  by 

matching  the  terms  of  Eq  (4.20),  namely  the  respective  terms  of  the 
equality 
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The  foilwing  pair  of  linear  equations  are  obtained 


U-XV  =  0 


(4.22) 


(U-XV)X^M-XN  =  M-XN  =  0 


(4.23) 


As  it  will  be  proved  below,  any  solution  of  Eq  (4,1)  must  satisfy  the 
above  pair  of  equations  (4.22)  and  (4.23).  In  other  words,  the  pair  of 
all  solutions  for  equations  (4.22)  and  (4.23)  trap  the  solutions  of  the 
Lyapunov  matrix  equation,  Eq  (4.1). 

Theorem  4-2.  If  N“^  exist,  and  X  is  a  common  solution  of 
equations  (4.22)  and  (4.23),  then  X  is  also  a  solution  to  Eq  (4.1). 

Proof.  Let  X  be  a  common  solution  to  Eqs  (4.22)  and  (4.23), 
then  using  the  identities  defined  in  Eqs  (4.16)  through  (4.19),  derive 
the  Lyapunov  matrix  equation  as  follows 

0  =  (XN-M)B 
=  XNB-MB 
=  XMB^UC-AM-CN 
=  X(BN-VC)-AM-CN+UC 
=  X(BN-VC)-AXN-CN^UC 
=  XBN-AXN-CN>(UC-XVC) 

=  (XBN-AXN-CN)>(U-XV)C 
=  (XB-AX-C)N 

■ 

Theorem  4-3.  Let  V-1  exist  and  let  X  be  a  common  solution  of 
Eqs  (4.22)  and  (4,23),  then  X  is  also  a  solution  of  Eq  (4,1). 

Proof.  Let  X  be  a  common  solution  to  Eqs  (4.22)  and  (4.23),  then 


using  the  identities  defined  in  Eqs  (4.16)  through  (4.19),  derive  the 
Lyapunov  matrix  equation  as  follows 


0  =  (U-XV)A 
=  UA-XVA 
=  XVA-(UA) 

=  XVA-(AU+CV) 

=  X(VA)-AU-CV 
=  XBV-AU-(CV) 

=  XBV-AU-(UA-AU) 
=  XBV-UA 
=  XBV-AU-CV 
=  -AU+XBV-CV 
=  -AXV-^XBV-CV 
=  (-AX+BX-C)V 


In  addition  to  the  two  Eqs  (4.22)  and  (4.23),  a  second  pair  of 
equations  can  be  used  to  trap  the  solution  to  the  Lyapunov  matrix 
equation.  The  common  solution  of  at  least  one  pair  will  define  the 
solution  to  the  Lyapunov  matrix  equation.  The  second  pair  of 
equations  are 


VX>N  =  0  (4.24) 

UX>M  =  0  (4.25) 


These  equations  are  derived  in  a  manner  similar  to  Eqs  (4.23)  and 
(4.24).  First,  the  equations  arise  from  the  results  of  theorem  4-1,  but 


with  two  terms  rearranged.  Now  the  similarity  transformation 
defined  for  the  matrices  of  Eq  (4.2) 


!‘"  and  ° 
[ij  bj  pj  5 


becomes  equivalently 


r  T  V’  1  [  i  1  i  T  _V 

I  J.  \  I  j.  1  •  I  1  i  i 


I J  [O  Bj  [0  I 


(4  26' 


Now  this  similarity  transformation  is  used  in  place  of  Eq  (4.20)  with 
fg(R)  substituted  for  fy^(R).  Now  similar  to  Eq  (4.20)  it  follows 


0  I 


rrl  Ft  ■'•.rl  Ft  ■'..’•1  Fit  Ft  _v 

"M  fp'iR)  M  w  1 


0  I  0  1  V  N  0  I 


I  X  U  M-UX 
0  I  V  M-VX 


(4.27) 


[F;+xy  (M-ux)+x(N-vx) 
I  V  N-VX 


I  -x'  ,  X^  Ip  (A)  0 

TplK.i  =  t' 

0  I  (i  I  *  fB®''' 


Here  fgfB)  is  equal  to  zero  which  provides  the  final  key  needed  to 
obtain  the  Eqs  (4,24)  and  (4.25).  For  a  more  complete  derivation  of 
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Eqs  (4.24)  and  (4.25)  see  (6:13-14).  Like  theorems  4-2  and  4-3,  Eqs 
(4.24)  and  (4.25)  tie  once  again  back  to  the  Lyapunov  matrix 
equation.  Theorems  4-4  and  4-5  show  this  with  the  given  conditions. 


Theorem  4-4.  If  exist,  and  X  is  a  common  solution  of 
equations  (4.24)  and  (4.25),  then  X  is  also  a  solution  to  Eq  (4.1). 

Proof.  Let  X  be  a  common  solution  to  Eqs  (4.24)  and  (4.25), 
then  using  the  identities  defined  in  Eqs  (4.16)  through  (4.19),  derive 
the  Lyapunov  matrix  equation  as  follows 

0  =  (VX^N)B 
=  VXB^(NB) 

=  VXB>BN-VC 
=  VXB-BV-VC 
=  VXB-VC-(BV)X 
=  VXB-VC-(VA)X 
=  V(XB-C-AX) 

■ 

Theorem  4-5.  If  U"^  exist,  and  X  is  a  common  solution  of 
equations  (4.24)  and  (4.25),  then  X  is  also  a  solution  to  Eq  (4.1). 

Proof.  Let  X  be  a  common  solution  to  Fqs  (4.24)  and  (4.25), 
then  using  the  identities  defined  in  Eqs  (4.16)  through  (4.19),  derive 


the  Lyapunov  matrix  equation  as  follows 


0  =  A(UX+N) 
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=  AUX^AM 
=  AUX^(AM) 

^  (AU)X+(UC+MB-CN) 

=  (UA-CV)X+UC+MB-CN 
=  UAX>UC-(CV)X^MB-CN 
=  UAX+UC-CN*MB-CN 
=  UAX+UC>(M)B 
=  UAX>UC>(-UX)B 
=  U(AX-XB*C) 


As  stated  earlier,  the  first  step  in  developing  a  solution  for  the 
Lyapunov  matrix  equation  was  deriving  an  intermediate  expression 
for  its  solution.  The  two  pairs  of  equations,  (4.22),  (4.23)  and  (4.24), 
(4.25),  comprises  this  intermediate  expression.  Now,  a  common 
solution  for  each  pair  must  be  obtained. 

Common  Solution.  Now  that  the  solution  of  the  Lyapunov 
matrix  equation  has  been  expressed  by  an  alternate  representation, 
the  final  step  is  to  find  the  solution  in  terms  of  this  representation. 
For  Eqs  (4.22)  through  (4.25),  this  means  finding  the  common  solution 
for  these  pairs  of  equations.  In  Odell  (31),  a  common  solution  for  n 
equations  is  derived  using  an  generalized  inverse.  Here  n  will 

be  fixed  at  two  and  an  Aj  2  will  be  substituted  in  place  of  the  Aj  2  3  4 
generalized  inverse  used  by  Odell. 

To  understand  why  an  A|  2  generalized  inverse  can  be  used  in 
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generalized  inverse  used  by  Odell,  examine  the 


place  the  Aj  2,3_4 

necessary  and  sufficient  condition  needed  for  a  common  solution  to 
exist.  Given  the  pair  of  equations 

AjX  =  Bj  (4.28) 

A2X  =  B2  (4.29) 

By  Odell  [31:272],  the  following  must  hold  for  Eqs  (4.28)  and  (4.29)  to 
have  a  common  solution 

^1^1(I,2,3,4)^I  " 

(A2‘A2AiAj)(A2-A2AiAj)(j,2,3,^)(B2-AiAjBj)  =  (4.31) 

(B2~AjA|B|) 


Examining  the  conditions  expressed  for  these  two  equations  reveals 
that  an  Aj  2  inverse  will  satisfy  both  conditions.  This  is  a  result  of 

the  partial  definition  of  an  Aj  2  generalized  inverse  which  was 
expressed  earlier  in  the  form 

AXA  =  A  (2.1) 

remembering  X  here  is  an  Aj  2  generalized  inverse.  For  Eq  (2.1)  to 

hold,  AX  must  equal  an  identity  matrix  of  appropriate  dimensions 
Applying  this  result  to  the  condition  stated  in  Eqs  (4.30)  and  (4.31) 
shows  that  the  Aj  2  satisfies  these  conditions  since  the  following 
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equations  also  hold 

^1^1(1, 2)®l  " 

M.32) 

(A2-A2AjAm2))(A2-A2AiAjQ2))(i2) 

(4.31) 

(B2"AjAjBj)  =  {B2“AjAjBj) 

since 

^1^1(1, 2)  "  ^ 

(A2-A2AiAj(i  2))(A2-A2AiAj(j  2))(i  2)  =  I 

by  the  definition  of  an  Aj  ^  generalized  inverse. 

Thus  using  an  Aj  2  generalized  inverse,  Odell’s  result  can  now 

be  stated  for  a  pair  of  equations  and  in  terms  of  Eqs  (4.24)  and  (4.25). 

Thus  given 

VX^N  =  0 

(4,24) 

UX>M  =  0 

(4,25) 

The  common  solutions  for  these  equation,  by  Odell  (31,272),  is  found 

from 

X  =  V,  2N^(I-Vi2V) 

(4.34) 

X  =  Vj  2NHI-Vj2V)[(U-UVi2V)j2(M-UVi2N)]  - 

(4.35) 

(I-Vj2V)[I-(U-UVi2V)j2  (U-UV12V)] 
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if  and  only  if  the  following  holds 


VVj  2N  =  N  (4.36) 

(U-Vj  2V)(U-Vi  2V)i  2(M-UVN)  =  (M-UVN)  (4.37) 

This  result  can  also  be  extended  for  Eqs  (4.23)  and  (4.24)  where 
the  unknown  is  on  the  left  hand  side.  In  this  case,  the  necessary  and 
suficient  conditions,  which  determine  whether  or  not  Eqs  (4.22)  and 
(4.23)  have  a  common  solution,  are  rearranged  to  account  for  the 
form  of  Eqs  (4.22)  and  (4.23).  Specifically,  the  second  characteristic  of 
an  Aj  2  generalized  inverse  is  utilized  to  obtain  the  new  expressions 

for  the  necessary  and  sufficient  conditions  and  also  the  form  of  the 
common  solutions  for  Eqs  (4.22)  and  (4.23).  This  characteristic  of  the 
Aj  2  generalized  inverse  was  stated  earilier  as 

XAX  =  X  (2.2) 

Now,  applying  Eq  (2.2)  to  the  conditions  needed  for  a  common 
solution  to  exist  for  Eqs  (4.23)  and  (4.24)  provides  a  similar  result  to 
that  developed  for  Eqs  (4.30)  and  (4.31).  This  result  can  be  stated  as 

V,  2VU  =  U  (4.38) 

(N-NV  Vj  2  )(N-NV  Vj  2  \2  (M-NUVj  2 )  =  (4.39) 

(M-NUV12) 

Now  if  Eqs  (4.38)  and  (4,39)  hold,  then  the  common  solution  for  Eqs 


(4.22)  and  (4.23)  can  be  obtained  from  the  following 

X  =  UVj  2*(I-VVi  2)  (4.40) 

X  =  UVj2*(1-UVi2)(M-NUVi2)(N-NVVj  2)1^2^ 

(I-Wi2)(I-(N-NVVj2)(N-NVVi  2)1^2) 


Probably  the  most  interesting  result  from  the  above  application 
of  a  generalized  inverse  to  the  Lyapunov  matrix  equation  is  the 
potential  savings  obtained  by  substituting  the  Aj  2  generalized  inverse 

over  the  Aj  2  3  4  generalized  inverse.  It  should  be  apparent  at  this 
point,  that  an  A|  2  is  easier  to  compute  than  an  Aj  2  3  4  (see  appendix 
C  for  a  quantitative  comparison  of  an  Aj  2  versus  an  Aj  2  3  4  )•  Thus  if 
an  Aj2  satifies  the  needs  of  an  application,  it  would  certainly  be 
more  desirable  than  an  Aj  2  3  4  for  computational  reasons.  This  is 

certainly  the  case  for  the  application  discussed  above  which  leads 
back  the  comments  made  in  the  introductory  remarks  to  this 
chapter  regarding  the  use  of  a  generalized  inverses.  Namely,  the  use 
of  a  lesser  generalized  inverse  was  probably  not  considered  for  solving 
for  the  common  solution  of  multiple  equations,  as  needed  above, 
primarily  due  to  the  lack  of  stable  and  direct  methods  for  computing 
such  an  inverse.  The  ST  method  provides  the  necessary  stable  and 
direct  computational  technique  needed  to  make  the  Aj  2  of  use  in 
this  application. 
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The  ST  method  provides  a  simple,  direct,  and  numerically  stable 
technique  for  computing  four  generalized  inverses  of  a  matrix.  It  also 
provides  a  unified  approach  for  computing  these  inverses  over  other 
techniques  which  normally  only  apply  for  computing  a  single 
generalized  inverse.  The  ST  method  also  highlights  the  computational 
differences  between  each  gereralized  inverse.  This  combined  with  the 
properties  of  a  respective  inverse  provides  insight  into  the  use  and 
cost  of  computing  a  generalized  inverse  for  a  given  application.  Thus, 
by  taking  into  account  the  cost  of  computing  an  inverse  and  its 
properties,  the  ST  method  may  help  reveal  where  a  lesser  generaliz¬ 
ed  inverse  will  satisfy  the  requirements  over  a  more  computational 
expensive  inverse.  This  was  certainly  the  case  for  finding  the 
solution  of  the  Lyapunov  matrix  equation.  In  summary,  the  ST 
method  is  a  simple,  direct,  and  unified  technique  for  computing  a 
generalized  inverse. 

There  are  two  areas  where  the  ST  method,  in  particular  its 
computer  implementation,  can  be  extended.  First,  further  study  is 
warranted  in  extending  the  ST  method  to  unbound  precision.  The 
next  logical  extention  would  be  to  address  matrices  with  single  or 
multiple  parameters.  Both  are  recommended  because  the  ST  method 
is  a  good  computational  framework  for  large  problems,  and  these 
potential  applications  deal  with  matrices  of  this  type  or  require 
these  capabilities. 


Appendix  A:  Example  1 


The  purpose  of  example  1  was  to  introduce  the  ST  method  in  a 
simple  and  direct  manner.  With  the  details  of  the  ST  method  now 
well  established,  example  3  extends  example  1  by  taking  a  more 
complete  look  at  the  steps  comprising  the  ST  method.  This  includes 
the  intermediate  results  of  the  row  and  column  operations  and  the 
use  of  Penrose's  four  equations  to  test  the  computed  generalized 
inverses.  The  matrix  A  in  this  example  was  taken  from  Noble 
(21:145, 146).  The  example  begins  with  the  given  matrix  A  of 
dimensions  6x4: 
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First  augment  A  with  an  identity  matrix  of  dimensions  6x6 
beside  A  and  an  identity  matrix  below  A  with  dimensions  4x4 
providing 
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Using  this  augmented  representation  of  A,  the  next  step  is  to 
reduce  the  original  matrix  to  a  square  identity  matrix  with 
dimensions  equal  to  the  rank  of  A.  The  reduction  process  uses  a 
modified  Guass- Jordan  algorithm  with  a  pivoting  scheme  which  will 
perform,  at  the  most,  three  pivots  to  locate  the  best  diagonal  element 
for  both  row  and  column  operations.  The  following  matrix  shows  the 
result  of  the  first  pivot  call  where  three  pivots  were  performed  to 
locate  the  best  diagonal  element  for  reducing  both  the  elements  below 
the  diagonal  and  the  elements  across  from  it: 
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With  the  best  diagonal  element  established,  the  first  column,  of  the 
now  shift  matrix  A,  is  eliminated  with  row  operations.  The  result  is 
shown  in  the  following  matrix; 
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The  pivot  process  continues,  but  now  the  pivot  is  performed  around 
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n:  I 

the  second  diagonal  element.  This  pivot  call  resulted  in  one  pivot 
which  is  shown  below; 
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Again,  the  next  step  is  to  reduce  the  elements  of  only  A  below  the 
diagonal  element  with  row  operations  which  produces 
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The  next  potential  diagonal  element  is  zero,  and  all  possible  pivot 
candidates  are  also  zeroes.  This  means  all  rows  below  the  last 
diagonal  element  have  been  reduced  and  no  further  row  operations 
are  needed.  This  also  denotes  that  the  T  submatrix  is  almost  defined 
since  it  is  a  function  of  the  row  operations.  The  only  remaining 
operations  which  may  impact  T  is  dividing  by  the  diagonal  element 
after  the  column  operations  have  been  completed.  This  will  be  more 
apparent  in  a  later  matrix. 

Since  the  pivots  performed  prior  to  each  set  of  row  operations 
were  designed  to  look  ahead  to  provide  a  good  computational  basis  for 
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the  column  operations,  no  f uther  pivoting  is  performed  during  the 
column  operations.  Instead,  the  column  operations  are  performed 
directly  on  the  remaining  non  zero  elements.  The  result  of  the  first 
set  of  column  operations  on  row  one  of  A  is  shown  below: 


3  000 

000 

000 

000 

000 

000 

1.000 

.000 

.000 

.000 

.000  -1 

000 

.  o’OsJ* 

.66? 

.000 

t.ooo 

.  333 
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.000 
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.000 

.000 

.000 
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.000 

.000 
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.000 

.  000 

.000 

.000 

.000 

.000 

1.000 
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.000 

.000 
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.  000 

.000 

.  000 

.  000 

.000 

.000 

.000 

1.000 

.000 

000 

1 . 000 

.000 

1 . 000 

000 

- .  333 

.  333 

The  next  set  of  column  operations  reduce  the  remaining  non  zero 
elements  of  row  two  across  from  the  second  diagonal  element.  This  is 
shown  in  the  following  matrix 


3 . 000 

.000 

.000 

.  000 

.000 

.000 
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.000 
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.000 

.000 

.000 
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.333 
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.000 
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.000 

.000 
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.000 
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.000 

.000 

.  000 
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.000 
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.000 

.000 
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.000 
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.667 
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.000 
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1 . 000 
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.  Osp'W 

With  only  diagonal  elements  remaing  in  A,  the  final  step  needed 
to  reduce  A  to  the  desired  identity  matrix  is  to  divide  through  by  the 
diagonal  elements.  With  all  row  and  column  operations  complete, 
the  submatrices  are  now  defined,  and  the  representation,  derived  in 
theorem  2-1  as  equivalent  to  the  augment  matrix  A,  is  obtained.  This 
representation  is  shown  in  the  next  equation  along  with  the  sub¬ 
matrices  o,  T,  M,  and  N  highlighted  with  the  partitioning  lines: 


•oo 


O 


1.000 

.000 

.  000 

1 . 000 

.  000 

.000 

.  000 

.  000 

.  000 

.000 

.000 

.  000 

.000 

1.000 

.000 

.000 

.  000 

.  000 

1 . 000 

.000 

33  .000 

000-1.000  -.333  .000 


1.000-1.000-1.000  .000 
.000  .000  1.000  1.000 


.000  1.000  .000 

.000  1.000  1.000 


.000  1.000 


.000  1.000 


.333  .667 

.000  1.000 


Using  the  above  matrix,  an  Aj  2  generalized  inverse  is  computed 

directly  from  the  product  of  the  S  and  T  submatrices.  This  product  is 
shown  in  the  following  matrix 


000 

- 1 . 000 

- .  333 

.  000 

.  000 

.000 

000 

.  000 

.  000 

.  000 

.  000 

.000 

000 

.  000 

.  000 

.  000 

.  000 

.000 

000 

.  000 

.  000 

.  000 

.000 

Applying  the  two  equations  an  Aj  2  generalized  inverse  satisfies 

provides  the  original  matrix  and  the  computed  inverse  as  follows  for 
each  respective  equation. 


1 . 000 

000 

1 . 000 

000 

1 .  CCIO 

1 . 000 

.  000 

-1 

000 

.  000 

- 1 . 000 

1 . 000 

000 

.  000 
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- 1 . 000 

“  J 

000 

1 . 000 

- 1 . 000 

.  000 

1 
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1 . 000 
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- 1 . 000 

-2 
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1,2^ 

^1,2  " 

^1,2 

000  -1 

.  000  - 

oc 

0 
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.  000 

.000 

.  000 
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.  000 

OtJO 

000 

000 

.  oc 

IJ 

000 
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■T'“'  0 

.  oc 

0 

000 
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'•vyv 


The  next  inverse  of  concern,  an  Aj  2  3  generalized  inverse,  is 

again  computed  from  the  product  of  the  S  and  T  submatrices,  but 
after  each  row  in  T  has  been  orthogonalized  to  every  row  of  M. 
Example  2,  in  chapter  3,  provides  more  details  regarding  the 
orthogonalization  technique  used.  Thus  an  Aj  2  3  generalized  inverse  is 

computed  from  the  following  framework  which  shows  the 
orthogonalized  T  and  M  submatrices  and  the  previous  S  submatrix. 


Ir|  0  I 
0  0  M 
SIN  0 


1 . 000 

.  000 

.  000 

.  000 

.056 

-.056 

.111 

-.111 

.056  -.056 

.  000 
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.056 
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.000 

.000 

.000 
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1 . 000 

.200  -.400 

000 

.  000 
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000 
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.000 

1.000 
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The  product  of  S  and  T  in  this  case  provides  the  following  matrix 
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This  Aj  2  3  satisfies  the  previous  two  equations  and  in  addition  the 
equation 


where  the  product  of  the  Aj  2  3  and  the  original  matrix  A  is  the 
symmetric  matrix  shown  below 
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To  obtain  an  Aj  2  4  generalized  inverse,  the  S  and  N  submatric¬ 
es  are  orthogonalized  like  the  T  an  M  submatrices  were 
orthogonalized  for  the  Aj  2  3  generalized  inverse  above.  The  product  of 
the  orthogonalized  S  submatrix  and  the  non  orthogonalized  T 
submatrix  form  the  Aj  2  4  generalized  inverse.  The  following  matrix 
shows  the  computation  framework  needed  to  compute  this  inverse; 
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The  product  of  S  and  T  in  this  case  produce 


000 

-647 

-  235 

000 

000 

000 

000 

.412 

.  059 

000 

.  000 

000 

000 

.235 

.  176 

.  000 

.  000 

.  000 

000 

.  059 

.  294 

.  000 

.  000 

.  000 

which  satifies  the  first  two  equations  and  the  equation 


The  product  of  Aj  2  4  and  A  again  produce  a  symmetric  matrix; 
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- .  059 
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235 
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.824 

The  final  generalized  inverse  an  Aj  2  3  4 ,  is  computed  from  the 

orthogonalized  version  of  both  S  and  T.  For  this  example,  the 
computational  framework  needed  to  compute  this  inverse  is 
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The  A  inverse  is 


-.  147 

- .  176 

.  029 

- .  029 
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.  127 
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.049 
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where  all  four  previous  equations  are  satisfied, 


DATE:  11/ 11/85 
VERSION;  10 

TITLE:  GENERALIZED  INVERSE  COMPUTATION  AND  VALIDATION 
LANGUAGE:  FORTRAN  77 

USE:  COMPUTE  AND  CKECKS  THE  COMPUTATION  OF  FOUR  GENERALIZED 
INVERSES. 


CONTENTS; 


MAIN  PROGRAM  -  INPUTS  MATRIX 

-  REDUCES  MATRIX  WITH  MODIFIED 
GUASS-JORDAN 

-  CONTROLS  TO  COMPUTE  GENERALIZED 
INVERSES 

I  PI  VOT  -  PERFORMS  A  THREE  STEP  PI  VOT 

MULT  -  MULTIPLIES  TWO  MATRICES  PLACING  THE 

RESULT  INTO  AN  OPTIONAL  THIRD  MATRX 
GROW  -  GRAM-SCHMIDT  COMPUTATION  FOR  ROWS 

GCOLUMN  -  GRAM-SCHMIDT  FOR  COLUMNS 

PRINTA  -  GENERAL  MATRIX  PRINT  ROUTINE 


FUNCTION: 

COMPUTES  THE  A,  2.  A  1^2,3,  ^1,2,4'  A|  2^3  4  GENERALIZED 

INVERSES  AND  VALIDATE  EACH  THE  FOUR  CONDITIONS  DEFINED 
BY  PENROSE. 


c 

C  DATE:  11/1 1/85 

C  VERSION:  10 

C 

C  NAME:  GINVERSE  (MAIN  PROGRAM) 

C  DESCRIPTION:  CONTROLS  THE  REDUCTION  OF  THE  INPUT  MATRIX  TO 
C  AN  IDENTITY  MATRIX  AND  ALSO  CALLS  ALL  NECESSARY 
C  SUBROUTINES  TO  COMPUTE  A  GIVEN  INVERSE 
C 

C  PASSED  VARIABLES:  A.ST,HOLD,R,C,ASIZE,RANK 
C  RETURNS:  NONE 
C  FILES  READ:  MATRIX  INPUT 

C  SUBROUTINES  CALLED:  IPIVOT,MULT,GROW,GCOLUMN,PRINTA 

C  CALLING  SUBROUTINE  NONE 
C 


PROGRAM  GINVERSE 


C  ARRAY  TYPE  AND  DIMENSIONS/  DECLARATIONS 


DOUBLE  PRECISION 
DIMENSION 

>  > 

p 

O 

DOUBLE  PRECISION 
DIMENSION 

ST 

ST(1:10.1:10) 

DOUBLE  PRECISION 
DIMENSION 

ORG 

0RG(1:10.1:10) 

DOUBLE  PRECISION 
DIMENSION 

HOLD 

H0LD(1:10.1:10) 

DOUBLE  PRECISION 
DIMENSION 

SR 

SR(1:10,1;10) 

C  VARIABLE  TYPE  DECLARATIONS 


INTEGER 

R 

INTEGER 

C 

INTEGER 

ASIZE 

INTEGER 

MINRC 

INTEGER 

RANK 

INTEGER 

1 

INTEGER 

J 

INTEGER 

K 

INTEGER 

TEST 

INTEGER 

OPTION 

INTEGER 

FORMATT 

INTEGER 

HEADER 

INTEGER 

SELECT 

INTEGER 

TEST 

INTEGER 

PIVOTRLT 

CHARACTER  FILENAM*10 

DOUBLE  PRECISION  MULTIPLIER 


PARAMETER(ZER0=5.E-I5) 

C  DATA  INITIALIZATION 


1000  FORMATOX;  SELECT  OUTPUT  OPTION’) 

1 00 1  FORMAT(  1 X;  OUTPUT  TO  SCREEN  -  SELECT  1  ’ ) 

1002  FORMATdX;  OUTPUT  TO  PRINTER  -  SELECT  2’) 

1003  FORMATdX) 

1004  FORMATdX;  SELECTION:') 

1005  F0RMAT(I2) 

1006  FORMATdX,’  SELECT  FORMAT  OPTION’) 

1007  FORMATdX,’  XX.XX  -SELECT  D 

1008  FORMATdX,’  XXX.XX  -SELECT  2’) 

1009  FORMATdX,'  XX.XXX  -SELECT  3’) 

1010  FORMATdX,’  XX.XXXX  -SELECT  4’) 

1011  FORMATdX,’  SELECT  INPUT  FILE’) 

1012  FORMATdX,’  CHOOSE  FROM  LIST  OR  SELECT  A  NEW  FILE:') 

1013  FORMATdX,’  1)  X.CASE  1’) 

1014  FORMATdX,’  2)  X.CASE  2’) 

1015  FORMATdX,’  3)  X.CASE  3’) 

1016  FORMATdX,’  4)  X.CASE  4’) 
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1017  FORMATdX;  5)  X.NOBLES  MATRIX’) 

1018  FORMATdX;  6)  NEW  FILE’) 

1019  FORMATdX,’  INPUT  FILE  NAME  IN  APPROPRIATE  FORMAT:  ’) 

1020  FORMAT(AIO) 

1021  FORMATCI*****  GENERALIZED  INVERSES  FOR  FILE:’, AlO,’ 

1 022  FORMATdX,’*****  ST  METHOD  LANGUAGE:  FORTRAN  77’) 

1023  F0RMAT(I2) 

WRITE(9,1000) 

WRITE(9,1001) 

WRITE(9,1002) 

WRITE(9,1003) 

WRITE(9,1004) 

READ  (9, 1005)  OPTION 

WRITE(9,1003) 

WR1TE(9,1006) 

WRITE(9,1007) 

WRITE(9,1008) 

WR1TE(9,1009) 

WRITE(9,1010) 

WRITE(9,1003) 

WRITE(9,1004) 

READ  (9, 1005)  FORMATT 

WRITE(9,1003) 

WRITE(9,1011) 

WRITE(9,1012) 

WRITE(9,1013) 

WRITE(9,I0I4) 

WRITE(9,1015) 

WRITE(9,1016) 

WRITE(9,1017) 

WRITE(9,1018) 

WRITE(9,1003) 

WRITE(9,1004) 

READ  (9, 1005)  SELECT 

IF(SELECT.EQ.1)THEN 

0PEN(  1  ,FILE="X.CASE  1”,STATUS=”0LD") 

FILENAM="X,CASE  1" 

END  IF 

IF  (SELECT.EQ.2)  THEN 


OPEN(  1  ,FILE="X.CASE  2",STATUS="0LD") 
FILENAM="X.CASE  2" 

END  IF 

1F(SELECT.EQ.3)THEN 

OPEN(  1  ,FILE="X.CASE  3“,STATUS="0LD") 
FILENAfi="X.CASE  3" 

END  IF 

IF  (SELECT.EQ.4)  THEN 

OPEN(  1  ,FILE="X.CA5E  5",STATU5="0LD“) 
FILENAM="X.CASE  5“ 

END  IF 

IF  (SELECT.EQ.5)  THEN 

OPEN(  1  ,FILE="X.NOBLES  MATRIX”,STATU5="0LD") 
FILENAM="X.NOBLES  MATRIX" 

END  IF 

IF  (SELECT.EQ.6)  THEN 
WRITE(9,1003) 

WRITE(9.1019) 

READ  (9,1020)  FILENAM 

0PEN(  1  ,F1LE=FILENAM.STATUS="0LD") 

END  IF 

WRITE(9,1003) 

WRITE(9,1003) 

IF(0PTI0N.EQ.1)THEN 
WRITE(9, 1 021)  FILENAM 
WRITE(9,I022) 

ELSE 

WR1TE(9, 1021)  FILENAM 
WRITE(9,1022) 

END  IF 

CREAD  SELECTED  FILE 

READ  (1,1023)  R 
READ  (1,1023)  C 
ASIZE  =  R+C 
DO  2001  1  =  1, ASIZE 
DO  2000  J=l, ASIZE 
READd,*)  A(I,J) 

ORG(l,J)  =  A(I,J) 

2000  CONTINUE 

2001  CONTINUE 


C  AUGMENT  A  WITH  THE  APPROPRIATE  IDENTITY  MATRICES 

DO  201 1  1  =  1, R 

D0  2010J=C+1,ASIZE 
TEST  =  J-C 
IF(TEST.EQ.I)THEN 
A(I,J)  =  1.0 
ELSE 

A(l,J)  =  0.0 
END  IF 

2010  CONTINUE 

2011  CONTINUE 

DO  2020  l=Rn,ASIZE 
DO  2021  J=1,C 
TEST  =  l-R 
IF  (TEST.EQ.J)  THEN 
A(I,J)  =  1.0 
ELSE 

A(l,J)  =  0.0 
END  IF 

2020  CONTINUE 

2021  CONTINUE 


C  DETERMINE  THE  MINIMUM  RANK 

IF  (R.LT.C)  THEN 
MINRC  =  R 
ELSE 

MINRC  =  C 
END  IF 

CALL  PRINT A(  1 , 1  ,ASIZE,ASIZE,  1 , OPTION, FORMATT.R.C, A) 
IF(0PTI0N.EQ.1)READ(9,«)XX 
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C  MATRIX  REDUCTION 


DO  3020  l=1,R-1 

PIVOTRLT=IPIVOT(R,C,MINRC,ASIZE.I,A) 
IF(PIV0TRLT.EQ.1)THEN 
D0  3010K=in,R 
MULTIPLIER  =  A(K,I)/ Ad, 1) 

A(K,l)  =  0.0 
DO  3000  J=I+I,ASIZE 

A(K,J)  =  A(K,J)-(A(I,J)’‘MULTIPLIER) 
3000  CONTINUE 

3010  CONTINUE 

ELSE 

GO  TO  3021 

END  IF 

3020  CONTINUE 

3021  CONTINUE 


3030  IF  ((ABS(A(l,l)).GT.ZERO),AND,(I.LT.C))  THEN 

DO  3050  J=IH,C 
MULTIPLIER  =  Ad, J)/A(l.l) 

Ad,J)  =  0.0 
DO  3040  K=RH,ASIZE 
A(K,J)  =  A(K,JHA(K,I)*MULTIPLIER) 
3040  CONTINUE 

3050  CONTINUE 
I  =h1 

GO  TO  3030 
END  IF 

RANK  =  0 
I  =1 
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3060  IF  ((ABS(A(l,l)).GT.ZERO).AND.(I.LE.MINRC))  THEN 
RANK  =  RANK+ 1 
MULTIPLIER  =  Ad. I) 

DO  3070  J=1.ASIZE 
A(I.J)  =  Ad.JVMULTIPLIER 
3070  CONTINUE 


GO  TO  3060 
END  IF 


CALL  MULT(R+ 1 , 1  ,C,  1  ,C*  1  .R,RANK,A,A,5T) 

CALL  PRINT A(  1 , 1  .C,R,3,OPTION,FORMATT.C,R,ST) 
IF  (OPTION.EQ.  DREAD  (9,DXX 


c 


CHECK  A|  2  USING 


1)  (AA^  2^^"^ 

CALL  MULT(  1 , 1  ,R,  1,1  .R,C,ORG,ST,HOLD) 

CALL  MULT(  1 , 1  ,R.  1 . 1  .C,R.H0LD,0R6,H0LD) 

CALL  PRINTA(  1 , 1  ,R.C.4.0PT10N.F0RMATT,R,C,H0LD) 
IF(0PTI0N.EQ.1)READ(9.‘‘)XX 

2)  2  =  ^1,2^ 

CALL  MULT(  1 . 1  ,C,  1 , 1  ,C,R,ST.0R6,H0L0) 

CALL  MULT(  1 . 1  .C.  1 , 1  ,R,C, HOLD, ST, HOLD) 

CALL  PRINT A(  1 , 1  ,C,R,5,0PTI0N,F0RMATT,C,R,H0LD) 
IF  (OPTION.EQ.  DREAD  (9,*)  XX 


C»«i*««C0MPUTE  A I  2,3  GENERALIZED  INVERSE 


4000 

4010 


DO  4010  l=),R 

DO  4000  J=CD,ASIZE 
SR(I,J-C)  =  A(I,J) 
CONTINUE 
CONTINUE 


CALL  GSROW(RANK.R,C,ASIZE,A) 

CALL  MULT(R+ 1 , 1  ,C,  I  ,C+ 1  ,R,RANK,A,A,ST) 

CALL  PRINT A(  1 , 1  ,C,R,6,0PTI0N,F0RMATT.C,R,ST) 
IF  (OPTION.EQ. DREAD  (9,*) XX 


CHECK  A 


1,2,3 


USING 


^^^1,2,3^*  "^^^1,2,3^ 


CALL  MULT(  1 , 1  ,R,  1 , 1  ,R,C,ORG,ST,HOLD) 

CALL  PRINT A(  1 , 1  ,R, R,7, OPT ION,FORMATT,R,R, HOLD) 
IF  (OPTION.EQ.  DREAD  (9,DXX 


'.t  TV'  r; » iv  zs  is  v  v  w  ii.^  m 


"v  'n.*w\*^.'  1^  TV  ' 


C«***«COMPUTE  A I  2,4  GENERALIZED  INVERSE 

CALL  GSCOL(RANK,R,C,ASIZE,A) 

CALL  MULT(R+ 1 , 1  .C,  U  ,R,RANK,A.SR,ST) 

CALL  PRINT A(  1 . 1 , C,R, 8, OPTION, FORMAT!, C,R, ST) 
IF(0PTI0N.EQ.1)READ(9.»‘)XX 


C  CHECK  A,  2,4  USING  (Aj  2,4  A)«  =  Aj  2,4  A) 

CALL  MULT( ) ,  1  ,C,  1 , 1  ,C,R,ST.ORG,HOLD) 

CALL  PRINT A(  1 , 1  ,C,C,9,0PTI0N,F0RMATT,C,C,H0LD) 
IF(0PTI0N.EQ.1)READ(9.*)XX 


C*»»**cOMPUTE  A I  2,3,4  GENERALIZED  INVERSE 


CALL  MULT(R+ 1 , 1  ,C,  1  ,C+ 1  ,R,RANK,A,A,ST) 

CALL  PRINT A(  1 , 1  ,C,R,  1 0,0PT10N,F0RMATT,C,R,ST) 
IF(0PTI0N.EQ.1)READ(9,»)XX 

C  CHECK  A|  2,3,4  USING  CONDITIONS  APPLIED  TO 

^1,2  ^1,2,3  ^1,2.4 


C 


CHECK  A,  2,3,4  WSING  (A  A,  2,3,4>*  =  A,  2,3,4 

CALL  MULT(  1 , 1  ,R,  1 , 1  ,R,C,ORG,ST,HOLD) 

CALL  PRINTA(1,1,R,R,1  l,OPTION,FORMATT,R,R,HOLD) 
IF  (OPTION.EQ. DREAD  (9,»‘)XX 


) 


C 


CHECK  A|  2,3,4  USING  (A|  2,3,4  A)*  =  A^  2,3,4 

CALL  MULT(  1 . 1  ,C, ) ,  I  ,C,R,ST,ORG,HOLD) 

CALL  PRINTA(  1 . 1  ,C,C.  1 2,0PTI0N,F0RMATT,C,C,H0LD) 
IF  (OPTION.EQ  DREAD  (9.*) XX 


% 


STOP 

END 


C  DATE:  11/11/85 

C  VERSION  10 

C 

C  NAME:  I  PIVOT 

C  DESCRIPTION:  FINDS  PIVOT  ELEMENTS  AND  SWITCHES  ROWS  AND/ 
C  OR  COLUMNS. 

C  PASSED  VARIABLES:  R,aMINRC,ASIZE,RANK,A 

C  RETURNS:  FUNCTION  VALUE 

C  FILES  READ:  NONE 

C  SUBROUTINES  CALLED:  NONE 

C  CALL  I NG  SUBROUT  I NE:  61 1 N  VERSE 

C 

FUNCTI0NIP1V0T(R,C,MINRC,ASIZE,CD,A) 


DOUBLE  PRECISION 

DIMENSION 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

DOUBLE  PRECISION 

DOUBLE  PRECISION 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 


A(1:10,l:10) 

R 

C 

ASIZE 

MINRC 

LARGEST 

EXCHANGE 

SWITCH 

SD 

SR 

SC 

I 

CD 

K 

COUNT 

ROW 

LASTSD 


PARAMETER(ZER0=5.E-15) 

SD  =  CD 
SR  =  CD 
SC  =  CD 
SWITCH  =  0 


mmmmmmm 


ROW  =  CD 
LASTSD  =  0 
COUNT  =  0 

1000  IF  ((SD.LE.MINRC).AND.(SWITCH.EQ.O.O))  THEN 

IF  (LASTSD.NE.SD)  ROW=SD 
IF  (LASTSD.EQ.SD).AND.(SC.EQ.CD)  60  TO  1021 
LASTSD  =  SD 

IF  (SC.EQ.CD)  LARGEST=ABS(A(SD,SD)) 

DO  1020K=SDn,R 

IF  (ABS(A(K,R0W)).LE.LAR6EST)  60  TO  1010 
LARGEST  =  ABS(A(K,ROW)) 

SR  =  K 
SC  =  ROW 

1010  CONTINUE 

1020  CONTINUE 

1021  COUNT  =  COUNT+ 1 

IF  ((LARGEST.GT.O.O).AND.(COUNT.GT.  1 ))  SWITCH=  1 
IF  ((SWITCH.EQ.  I  ).AND.(C0UNT.6T,  1 ))  GO  TO  1050 

DO  1040  l=SDH,C 

IF  (ABS(A(SD,I)).LE.LARGEST)  GO  TO  1030 
LARGEST*  ABS(A(SD,  I )) 

SC  =  I 
SR  =  SD 
ROW*  I 

1030  CONTINUE 

1040  CONTINUE 

1050  CONTINUE 

IF  ((SWITCH.EQ.0).AND.(SC.EQ.CD),AND.(5R.EQ.CD))  SD=SDH 
GOTO  1000 
END  IF 


IF  ((SC.EQ.CD).AND.(SR.EQ.CD).AND.(LARGEST.GT.O.O))  THEN 


IF  (SC.NE.CD)  THEN 
DO  1060  l  =  l,ASIZE 
EXCHANGE  *  Ad, SC) 
A(I,SC)  =  A(I,CD) 
Ad, CD)  =  EXCHANGE 
CONTINUE 
END  IF 


IF(SR.NE.CD)THEN 
DO  1070  l=l,ASIZE 
EXCHANGE  *  A(SR,1) 
A(SR,I)  =  A(CD,I) 
A(CD,I)  =  EXCHANGE 
CONTINUE 
END  IF 


IF  (ABS(A(CD,CD)).GT.ZERO)  THEN 
IPIV0T=1 
ELSE 
[P(V0T=0 
END  IF 


RETURN 


mT  ^  , 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DATE.  11/11/85 
VERSION:  10 


NAME:  PRINTS 

DESCRI PTION:  PRINTS  HEADER  VARIABLE  HEADER,  FORMAT  AND 
MATRIX  DATA  BASED  ON  V ARI ABLES  PASSED 
PASSED  VARIABLES:  STR,  STCER,EC, HEADER, FORMATT,R,C.A 
RETURNS:  NONE 
FILES  READ:  NONE 
SUBROUTINES  CALLED.  NONE 
CALLING  SUBROUTINE:  GINVERSE 


SUBROUTINE  PRINT A(STR,STC,ER,EC.HEADER,OPTION,FORMATT,R,C,A) 


DOUBLE  PRECISION  A 
DIMENSION  Ad:  10, 1:10) 


INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 


R 

C 

STR 

STC 

ER 

EC 

HEADER 

OPTION 

FORMATT 


PARAMETER(ZER0=5.E-I5) 

1000  FORMAT(  IX, 'Initial  matrix  with  augmented  Identity  matrices') 

1001  FORMAT(  IX, 'Reduced  matrix  with  S,  T,  M,  and  N  submatrices') 

1002  F0RMAT(1X,'A1,2  generalized  inverse') 

1003  FORMAT(lX,'Check  Al,2  against  the  following  conditions:') 

1004  FORMAT(  IX,'  (A  A1 ,2  A)  =  A') 

1005  FORMATdX,'  (Al,2  A  Al,2)  =  Al,2') 

1006  F0RMAT(1X,'A1,2,3  generalized  inverse') 

1007  FORMAT(lX,'Check  A  1,2,3  against  (A  Al,2.3)*  =  (A  Al,2,3)') 
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1008  FORMAK IX, ‘A  1,2, 4  generalized  inverse*) 

1009  FORMAK lX,*Check  A  1,2,4 against  (Al,2,4  A)*  =  (A1,2,4  A)*) 

1010  FORMAK  1 X,’ A 1, 2,3,4  generalized  inverse*) 

1011  FORMAK  lX,*Check  A1 ,2,3,4  with  (A  A1 ,2,3,4)*  =  (A  A1 ,2,3,4)*) 

1012  FORMAK  1  X,*Check  A 1 ,2,3,4  with  ( A1 ,2,3,4  A)*  =  ( A 1 ,2,3,4  A)*) 

1013  FORMAK  IX) 

1014  FORMAK  12(F5.2,  IX)) 

1015  FORMAK  12(F6.2,  IX)) 

1016  FORMAK  12(F6.3,  IX)) 

1017  FORMAK  12(F7.4,  IX)) 

1F(0PTI0N.NE.1)G0T0  2000 

IF(HEADER.EQ.l)  WRITE(9,1000) 

IF  (HEADER.EQ.2)  WRITE(9,1001 ) 

IF  (HEADER.EQ.3)  WRITE(9,1002) 

IF  (HEADER.EQ.4)  WRITE(9,1003) 

IF  (HEADER.EQ.4)  WRITE(9,1013) 

IF  (HEADER.EQ.4)  WRITE(9,1004) 

IF  (HEADER.EQ.5)  WRITE(9,1005) 

IF  (HEADER.EQ.6)  WRITE(9,1006) 

IF  (HEA0ER.EQ.7)  WRITE(9,1007) 

IF  (HEADER.EQ.8)  WRITE(9,1008) 

IF  (HEADER.EQ.9)  WRITE(9,1009) 

IF  (HEADER.EQ.  1 0)  WRITE(9, 1010) 

IF  (HEADER.EQ.  1 1 )  WRITE(9, 1011) 

IF  (HEADER.EQ.  12)  WRITE(9,1012) 


GOTO  2010 

IF  (HEADER.EQ. 
IF  (HEADER.EQ. 
IF  (HEADER.EQ. 
IF  (HEADER.EQ. 
IF  (HEADER.EQ. 
IF  (HEADER.EQ. 
IF  (HEADER.EQ. 
IF  (HEADER.EQ. 
IF  (HEADER.EQ. 
IF  (HEADER.EQ. 
IF  (HEADER.EQ. 
IF  (HEADER.EQ. 


1)  WRITE(10,1000) 

2)  WRITE(10,1001) 

3)  WRITE(10,1002) 

4)  WRITE(10,1003) 
4)  WRITE(10,1013) 

4)  WRITE(10,1004) 

5)  WRITE(10,1005) 

6)  WRITE(10,1006) 

7)  WRITE(10,1007) 

8)  WR1TE(10,1008) 

9)  WRITE(10,1009) 

10)  WRITE(10,1010) 


no 


IF  (HEADER.EQ.  11 )  WRITE(  1 0. 1 0 11 ) 
IF  (HEADER.EQ.12)  WRITE(10.1012) 


2010  CONTINUE 


IF(0PTI0N.EQ.1)THEN 

WRITE(9,1013) 

ELSE 

WRITE(10,1013) 
END  IF 


2020 


DO  2060  l=STR,ER 
DO  2020  J=STC,EC 

IF  (ABS(A(l,J)).LT.ZERO)  A(I.J)  =  0.0 
CONTINUE 


2030 


IF(0PTI0N.NE.1)G0  TO  2030 
IF  (FORMATT.EQ.  1 )  WRITE  (9, 1 0 1 4)  (A(I,J),J=STC,EC) 
IF  (F0RMATT.EQ.2)  WRITE  (9,1015)  (A(I,J),J=STC,EC) 
IF  (F0RMATT.EQ.3)  WRITE  (9,1016)  (A(I,J),J=STC,EC) 
IF  (F0RMATT.EQ.4)  WRITE  (9,1017)  (A(I,J),J*STC,EC) 

GO  TO  2040 

CONTINUE 

IF  (F0RMATT.EQ.1)  WRITE(10,I014)  (A(I,J),J=STC,EC) 
IF  (F0RMATT,EQ.2)  WRITE(10,1015)  (A(I,J),J=5TC,EC) 
IF  (FORMATT.EQ, 3)  WRITE(10,1016)  (A(I,J),J=STC,EC) 
IF  (F0RMATT.EQ.4)  WRITE(10,1017)  (A(I,J),J=STC,EC) 


2040  CONTINUE 


2060  CONTINUE 


DATE;  11/11/85 
VERSION:  10 


NAME:  MULT 

DESCRIPTION:  MULTIPLIES  ANY  TWO  MATRICES  OR  SUBMATRICES 
AND  PLACES  THE  RESULT  IN  AN  OPTIONAL  THIRD 
MATRIX  OR  SUBMATRIX. 

PASSED  VARIABLES;  XR,  XC.  Dl,  YR,  YC,  D2,  CONFORMITY,  Ml 

M2.  AND  M3  * 

RETURNS;  NONE  * 

FILES  READ:  NONE 

SUBROUTINES  CALLED:  NONE  * 

CALLING  SUBROUTINE  GINVERSE 


SUBROUTINE  MULT(XR,XC,D  1  ,YR.YC,D2.CONFORMITY,M1  ,M2,M3) 

DOUBLE  PRECISION 

Ml, M2, M3 

DIMENSION 

MI(1:10,I:I0) 

DIMENSION 

M2(1:I0,1;10) 

DIMENSION 

M3(1:10.l:10) 

INTEGER 

XR,XC,D  1  ,YR,YC,D2,CONFORMITY 

DOUBLE  PRECISION 

ROW 

DIMENSION 

ROW(l:IO) 

INTEGER 

DXR 

INTEGER 

DXC 

INTEGER 

DYR 

INTEGER 

DYC 

INTEGER 

IXR 

INTEGER 

JYC 

INTEGER 

KXC 

INTEGER 

KYR 

INTEGER 

1 

INTEGER 

J 

INTEGER 

K 

DXR*XR-1 

DXC  =  XC-1 

DYR»YR-1 

DYC  =  YC-1 


1010 

1020 


1030 


DO  1040  1=1, D1 
IXR  =  l+DXR 

DO  1020J=1,D2 
JYC  =  J+DYC 
R0W(J)=0.0 

DO  1010  K=1, CONFORMITY 
KXC  =  DXC+K 
KYR  =  DYR>K 

ROW(J)  =  R0W(Jh(M1(IXR.KXC)*M2(KYR  JYO) 
CONTINUE 
CONTINUE 
DO  1030J=1,D2 
M3(I,J)  =  ROW(J) 

CONTINUE 


1040  CONTINUE 


RETURN 

END 


m 


DATE;  11/11/85 
VERSION:  10 


NAME:  GSROW 

DESCRIPTION:  ORTHOGONAL  I ZES  THE  ROWS  OF  THE  SUBMATRICES 
T  ANDM 

PASSED  VARIABLES;  RANK,  R,  C,  ASIZE,  A 
RETURNS:  NONE 
FILES  READ:  NONE 
SUBROUTINES  CALLED;  NONE 
CALLING  SUBROUTINE;  GINVERSE 

SUBROUT  I NE  GSROW(RANK,R,C,  ASI ZE,  A) 

DOUBLE  PRECISION  A 
DIMENSION  Ad:  10, 1:10) 

INTEGER  RANK,R,C,  ASIZE 

DOUBLE  PRECISION  D0T1 
DOUBLE  PRECISION  D0T2 
INTEGER  I 

INTEGER  J 

INTEGER  K 


PARAMETER(ZER0=5.E-15) 


DO  1040  l=R,RANKd,-l 
DOT2=0.0 

DO  IOOOK=Cd,ASlZE 

D0T2  =  D0T2*(A(I,K)**2) 
CONTINUE 
DO  1030  J=l-1, 1,-1 
DOT  1=0.0 

DO  1 01  OK=Cd, ASIZE 
DOTl  =DOTU(A(l,K)*A(J,K)) 
CONTINUE 


IF  ((ABS(D0T1).GT.ZER0).AND.(ABS(D0T2).GT,ZER0))  THEN 
DOTl  =D0T1/D0T2 
DO  1020  K=C+1.ASIZE 
A(J,K)  =  A(J.KHD0T1*A(I,K)) 

1020  CONTINUE 

END  IF 

1030  CONTINUE 
1040  CONTINUE 


RETURN 


c 

C  DATE:  11/11/85 

C  VERSION  1.0 

C  * 

C  NAME;  6SC0LUMN 

C  DESCRIPTION:  0RTH060NALIZE  THE  COLUMNS  OF  THE  SUBMATRICES  * 

C  S  AND  N 

C  PASSED  VARIABLES;  RANK,  R.  C,  ASIZE,  A 

C  RETURNS:  NONE  * 

C  FILES  READ;  NONE 

C  SUBROUTINES  CALLED:  NONE 

C  CALLING  SUBROUTINE:  NONE 

C  * 
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Appendix  C:  ^kaph.  oi  ^  Ai^2  ^^^'sus  an  Ai^2,3,4 


v.v 


Starting  with  the  A^  ^  continuing  through  the  A^  2  3  4 . 

each  successive  inverse  satisfies  either  more  or  a  different  combina¬ 
tion  of  Penrose's  equations.  This  relationship  between  each  generaliz¬ 
ed  inverse  is  also  reflected  in  the  computer  time  needed  to  compute 
an  inverse.  To  meet  the  additional  constraints,  successive  inverses 
require  the  same  or  more  computer  time.  This  relationship  was 
originally  mentioned  in  chapter  3  and  its  practical  aspects  mentioned 
again  in  chapter  A  where  an  Aj  2  used  in  place  of  an  Aj  2^3^4 
generalized  inverse.  The  following  graph  emphasizes  this  realtionship 
by  highlighting  the  differences  in  computing  time  between  an  Aj  2 

and  an  Aj  2  3  4  over  a  range  of  square  matrices. 
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In  this  work,  four  generalized  inverses  are  derived  and 
computed  in  a  systematic  manner  from  this  representation.  These 
results  are  strongly  tied  to  the  solution  of  the  matrix  equation 
Ax=b  where  the  general  solution  is  given  in  terms  of  this  new 
representation  and  computational  technique.  The  computational 
technique  is  presented  with  an  example  and  also  as  an  algorithm. 
Included  with  the  algorithm  is  an  analysis  of  its  computer 
implementation.  The  use  of  one  generalized  inverse  is  used  to  find 
the  solution  of  a  Lyapunov  matrix  equation. 
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